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Summary. The two-dimensional convective motion generated by buoyancy forces 
on the fluid in a long rectangle, of which the two long sides are vertical boundaries held 
at different temperatures, is considered with a view to the determination of the rate of 
transfer of heat between the two vertical boundaries. The governing equations are set 
up; they reveal that the flow is determined uniquely by the Prandtl number a, the 
Rayleigh number A = g(T,; — T,)d*/(Toxv), and the ratio of the sides of the rectangle 
l/d. In the case of cavities used for thermal insulation of buildings, which is kept specially 
in mind throughout the paper, A is usually about 1000 d* (where d is in centimeters), 
and //d takes values between about 5 and 200. 

The essence of the problem is to determine which of several different flow regimes 
occurs at any given values of A and //d. It appears that with the above practical values 
the flow is not decisively of one single kind, and the discussion of the heat transfer for 
each of several ranges of values of A and I/d is necessary. Sections 4, 5 and 6 are con- 
cerned with laminar flow regimes characterized by very small values of A, large values 
of //d, and large values of A, respectively. In Section 7 approximate criteria for these 
flows to be stable are established, and in Section 8 the expressions for the heat transfer 
when the flow is turbulent are considered briefly. The unified picture provided by all 
these different results is considered in Section 9. 

A comparison of the theoretical predictions about the heat transfer with the limited 
experimental data, mostly obtained by Mull and Reiher (1930), is made. Theory and 
experiment agree in suggesting that, under practical conditions, the effect of convection 
is negligible for d < 1 cm, and that the heat transfer per unit area of vertical boundary 
decreases as d increases, provided d < 2.5 cm, and remains approximately constant 
(at a value proportional to /~'”*) for further increase of d. 

1. The background to the problem. The purpose of this paper is to determine the 
rate at which heat is transferred across the air space between two plane parallel vertical 
boundaries which are held at different temperatures. The air space is closed by horizontal 
boundaries distance | apart (l being large compared with the distance d between the 
vertical walls, in general), as sketched in Fig. 1. In the remaining direction, at right 
angles to the plane of the sketch, the air space is regarded as extending to infinity. All 
boundary conditions will be assumed uniform in this latter direction, and the convective 
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motion generated by buoyancy forces can then be assumed to be two-dimensional, lying 
in the plane of Fig. 1. Only the heat transferred from one vertical boundary to the other 
by conduction and convection will be considered. The radiative transfer is not negligible, 
but it takes place independently of the conductive and convective transfer and does 
not depend significantly on the shape and size of the air space; once the nature of the 
surface of the boundaries and their temperature difference have been specified, the 
radiative transfer can be calculated with reasonable accuracy,from known laws. 
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The determination of the transfer of heat across an air space of given size and shape 
is a problem which arises frequently in connection with the thermal insulation of build- 
ings. This is the context in which the problem first came to my notice, and it will be kept 
in mind in the analysis that follows. It has long been appreciated by those concerned 
with the thermal insulation of buildings that a narrow gap or cavity in the interior of a 
wall can impede the flow of heat considerably without adding appreciably to the cost of 
construction of the wall. For instance the rate of transfer heat through a brick wall 9 
inches thick is about 40% greater than that (with the effect of radiation included) 
through a composite wall consisting of two 43” brick leaves with an unventilated 
cavity 2 inches wide between them, given the same temperature difference between the 
two outer brick faces in each case*. The use of cavity walls of this kind has been standard 
practice for some years. Double windows consisting of two panes of glass with an air 
space between them are also well-known in principle as a means of improving thermal 
insulation (although it does not seem to be so widely appreciated in Great Britian that 
their use may also be economically desirable, at any rate for rooms which are kept at a 
comfortable temperature throughout most of the year). 

Now it is clear that if the air space between the two vertical boundaries is very 
narrow, very little convective motion can occur and the heat transfer will be due mostly 
to conduction. If the air gap is widened, the transfer that would result from conduction 
alone will decrease but the convective transfer will increase. Thus there exists the 
possibility that the transfer is a minimum for an air space of certain width. Inasmuch as 
the convective motion is mostly in the vertical direction, it is likely that the height of 
the air space also has an influence on the heat transfer. In the design of a building in- 


*An excellent account of the principles and practice of thermal insulation of buildings (with English 
conditions chiefly in mind) is given by N. 8S. Billington, in ‘“Thermal] Properties of Buildings,’’ Cleaver- 
Hume Press Ltd. 1952. 
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volving the use of an air space as a thermal barrier the gap d, and to a lesser extent the 
height / will be disposable quantities, unlike the temperature difference, so that the 
practical value of the investigation will lie in the information it provides about the 
dependence of the heat transfer on d and l. This information is of special importance for 
a consideration of the use of double windows, since first, the thermal resistance provided 
by the air space is here more than half the total resistance of the window, and second, 
the use of double windows will usually be economically justified by a small margin only, 


\lthough no anlytical discussion of the determination of the heat transfer seems to 
have been published, a number of measurements of the heat transfer have been made; 
the relevant parts of the data will be described later. So far as I can gather from the 


literature of heating engineers, it seems to be generally accepted that an optimum size 
of air gap exists and that it is close to 3/4 inch, at any rate for double windows. The fact 
that the heat transfer may depend on the height | seems not to be well appreciated, and 
the possibility that the optimum value of d itself depends on / has not always been taken 
into account in the design of the experiments. No observations of the velocity or tempera- 


ture distributions within the air space seem to have been published. 

2. The equations governing the problem. Let 7, and 7, be the absolute tempera- 
ures of the two vertical boundaries. The main assumption underlying the equations 
to be used is that the temperature difference 7, — 7’) is a small fraction of the absolute 
temperature 7, and that the variation of temperature in the fluid can be neglected for 
all purposes other than the determination of the buoyancy force. The self-consistency 
of equations based on this assumption for problems of free convection has been explained 


elsewhere (Goldstein, 1938). We shall see in the next section that (7, — 1T,)/To> is in 
fact small in eases of heat transfer in buildings. 
Now it is clear that the speeds involved in problems of free convection of this kind 


are far below the speed of sound, and that the pressure differences produced by inertia 
forces are a minute fraction of the absolute pressure. The pressure differences produced 


by gravity are also very small compared with the absolute pressure (being controlled 
by the length scale of the space occupied by fluid), and variations of the fluid density 
p will be determined wholly by variations of the temperature 7’. The equation of state 
for gases then has the form 
> = i ans Zo 9 
oe ees 2.1) 
Po To 


this and subsequent equations can be made applicable to liquids also by replacing the 
factor 7 on the right hand side by a coefficient of expansion, since the magnitude of 
the temperature does not occur in the equations in any other connexion. 

Che equation expressing conservation of mass is the same, with the assumptions 
described, as that for an incompressible uniform fluid. Thus if coordinates x, y are 
chosen as in Fig. 1, and the corresponding velocity components are u, v, we have 


Ou Ov — 29 
Ox Oy 


The only cause of change of the temperature of a moving element of fluid is heat 
conduction in the fluid, since neglect of temperature changes due to compression of the 
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fluid and due to heat of viscous dissipation is consistent with the assumptions already 
made. Hence, on introducing 


i igs Te 
9 = —_——* 
T; ~~ To 
as a more useful variable than 7, we have 
28 204 58. (28 4 20) - 
at +4 T°” "\act * ay? (2.3) 


where « is the thermal diffusivity (= k/poc, for gases, k/poc, for liquids, where k is the 
thermal conductivity of the fluid, and c,, c, are the specific heats); x is effectively uniform 
in view of the smallness of (7, — 7T>)/To. 

The force equations, taking account of the effect of buoyancy, are 


du, du, ou_ _ 1 dp (2 = t») (a *) 9 

at reat dy sp +9 Po T+ aa + oy’ /’ (24) 
wy, _ _ il op (2, #) 9 5 
at i dx vs dys po dy \Oa” + ay"/’ (2.5) 


where p is the pressure in the fluid. In these equations the density and kinematic viscosity 
v have been taken as uniform (and equal to p, and yuo/p. respectively, where yu is the 
fluid viscosity at temperature T,), again in view of the smallness of (7, — 7>)/T> 

The boundary conditions that will be imposed on the variables u, v and 6 at the vertical 


boundaries are 


u=v=0 6 = 0, at y = 0, 


’ 


u=v=0, 6 = 1, at y = d, 


the assumption being that these boundaries are made of material of high conductivity 
not in the form of thin sheets. On the two sides of the boundary given by x = 0 and 
xz = Il, we shall require u = v = O, together with a reasonable condition on @. In cases 
in which the boundary is of the same material on all four sides, the only consistent simple 
assumption is that 6 = y/d on the horizontal boundaries. However in cases in which the 
vertical and horizontal boundaries are not both appreciably better conductors than the 
fluid (as for instance where air is enclosed by two glass panes set in a wooden frame; 
Keines = 30K,;,, and kyoog = 6k,;,), this assumption may not be accurate, and it may be 
useful to consider the other extreme case of insulating horizontal boundaries, for which 
the boundary condition is 0@/8x = 0. The length of horizontal boundary is small com- 
pared with the length of vertical boundary, and it seems unlikely that the exact form 
of the boundary condition at x = 0 and x = I has much effect on the heat transferred 
through the vertical boundaries. 

These are the general equations and boundary conditions that will govern all types 
of flow in the cavity. Since the time scale of changes in temperature inside or outside 
buildings is large compared with time scales related to the convective motion, steady 
motions (or steady mean motions in the case of turbulent flow) are of interest herein 
and we therefore put 8/dt = 0. Further simplification is obtained by writing the equa- 
tions in dimensionless form. We take d as the unit of length for the coordinates z and y 
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(without changing the notation), and define a dimensionless stream-function ¥, with 
the aid of (2.2), by the relations 


gh 
u = xd ay ? v = —«xd ro (2.6) 
The heat equation (2.3) can then be written as 
addy daddy _ A(6, W) = 76. (2.7) 


dx dy dydx Az, y) 


When the pressure p is eliminated from (2.4) and (2.5), and p is eliminated with the use 
of (2.1), we are left with the second governing equation for y and @: 


1d@, W _ , 98 2 
22,9) Aa + V™ as 
where w = — 7’ is the (dimensionless) vorticity, ¢ is the Prandtl number v/«x, and 
— 3 
A = i= ted (2.9) 


T Kv 


is the Rayleigh number*. The boundary conditions are now 
= Ps _ 
=? 6=0, at y = 0, 
oy _ 
oy 


y= =0, @=y or ad at z=0 and t= 


Qs 


The form of the governing equations (2.7) and (2.8) and the above boundary con- 
ditions shows that the dimensionless parameters whose values are sufficient to determine 
uniquely the distributions of y and @ are o, A and l/d. 

The quantity to be determined from the analysis is the rate at which heat is trans- 
ferred by conduction through either vertical boundary. If this rate is Q heat units per 
second per unit depth of boundary (in the z-direction), a suitable dimensionless quantity 
describing the heat transfer is the Nusselt number 


Q * (20) 

a =) 2.10) 
KT, — To) Jo \dY/y-0 

N has the value //d when there is no convection, and in the general case we anticipate 

that 





N = N(o, A, U/d). (2.11) 


*The Grashof number G = A/c is frequently used by authors in the analysis of problems of free 
convection, but the number A enters directly at least as often as G. Moreover, as Professor H. B. Squire 
has pointed out to me, there is a much stronger case for honouring the name of Rayleigh in this field. 
The symbol A is used here to avoid the introduction of either R, which usually denotes Reynolds number, 
or the double letter Ra. 
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A quantity which is perhaps more useful for practical design of buildings is the heat 
transfer coefficient or thermal conductance of the cavity, defined by 


Q a 
C a tele Et 2.12) 
C is the rate of heat transfer through the vertical boundary, per unit area, per unit 
temperature difference between the boundaries. 
It can be seen immediately that if the temperature in the cavity were distributed 


as though the air were at rest, i.e. 0 = y, the term representing the effect of buoyancy in 
Eq. (2.8) is finite, and it is not possible for both of the other two terms to vanish. Zero 


velocity everywhere is thus not a possible solution of the equations, and any temperature 


difference across the cavity, however small, will produce some steady convective motion. 
The problem is quite different from that in which the temperature gradient is across the 
two horizontal sides of a cavity, when steady convective motion can occur only if the 
Rayleigh number exceeds a certain critical value. A plausible inference from the existence 


at all Rayleigh numbers is that the dependence of y and @ on A 


of a convective motior 
is smooth and that it is possible to expand y and @ as power series in A for sufficiently 
his provides a potential method of solution of the problem which 


small values of A. Th 
will be taken up in section 4, although it will be found that the first few terms of this 


series provide a valid approximation only when the convective flux of heat is small 


r 


ransfer by conduction, and its usefulness is very limited. 


compared with the t 
[t will be useful to have some idea of the values of A 


3. The practical conditions. 
appropriate to the practical problem of heat transfer through the walls of buildings. 
For air at 10°C and atmospheric pressure we have 


= ().14 cm’ sec , x = 0.19 em’ sec, ¢ = 0.73. 


The difference between the temperature of the air inside the building and that outside 
PS 
and an average value during the winter 


will not often exceed 20°C in Great Britain 
15°C for an indoor temperature of 18°C (= 65°F)*. Of this 


months would be about 
C.. (where C,, is the air-to-air conductance 


air-to-air temperature difference, a fraction C 
—including the effect of radiation—of the whole composite wall and C, is the boundary- 
to-boundary conductance of the cavity) represents the difference between the tempera- 
tures of the boundaries of the cavity. For double windows C,,/C, is roughly 0.6, and for 
a wall composed of two 44” brick leaves with a cavity C,,/C, is roughly 0.4. Taking 


the value C,,/C. = 0.5 we have 7.5°C as a representative value of 7, — 7%, and the 
corresponding value of A for d = 1 em is (very conveniently) 
980 X 7.5 ’ 
A, = - = 1000. (3.1) 


276 X 0.19 X 0.14 


The values of d used in practice range from about 1 cm up to about 8 em and the corre- 
sponding values of A for the conditions described by (3.1) (which will be referred to as 
the “standard” conditions) are given by 


A = 1000d’. (3.2) 


where d is expressed in centimetres. 


*The temperature difference appropriate to conditions in North Europe and North America might 


well be twice as large. 
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For double windows 1 will normally lie between 25 and 200 cm, and values of l/d 
ranging from about 5 to 200 should be considered in the analysis. In the case of a cavity 
in a brick wall, J will normally be at least 250 cm and may be much larger; standard 
building practice is to make the cavity 5 to 7 em wide, giving a range of values of I/d 
which is included within that just mentioned. 

4. The solution for very small values of A. We consider here a method of solution 
by means of an expansion of the quantities @ and y in power series in A. As already 
intimated the method is useful only at values of A too small for appreciable convection 
to occur, but it is necessary to give a few details of the method in order to be able to 
estimate its range of validity. 

The expansions to be employed are 


A(x, y) wat + A6,(x, y) + A’6,(x, y) + °° ’ (4.1) 
V(x, y) sd Ay,(z, y) + A* p(x, y) + °° ; (4.2) 


and the boundary conditions on the coefficients are 


_ In _ a : = , = 
y, = * 6, = 0 at, y=0 and y= 1, 





«= 2% 0, 6, = 0 or we at z=0@ and rai. 


These series can be substituted into equations (2.7) and (2.8), whence by equating 
coefficients of like powers of A the following set of equations is obtained; 
dy 
‘A =1 = —-— 4. 
V V1 ’ V 6, Ox ’ ( 3) 
_ 96 , 1HV*h , Ws) 2p — — 942 , 96, Wi) 
— Oy a Ax, y) ’ van ~ | * A(z, y) ’ ated 


ee 


The distribution of y, is thus identical with the distribution of displacement of an elastic 
plane rectangular plate clamped at the edges and subjected to a (small) uniform trans- 
verse pressure. An analytic solution to this latter problem is not known, so that it is 
necessary to determine y, numerically. When that has been done, the determination of 
@, requires a numerical integration (which can be carried out explicitly), and y. , 4, , etc., 
can then be determined by similar procedures. 

It is clear from the form of Eqs. (4.3), (4.4), ete. that y, and @, have the 
same symmetry properties as (x — h/2d)"~' (y — 1/2)""" and (x — h/2d)" (y — 1/2)""* 
respectively. Hence the horizontal flux of heat is given by 


l/d a © 
[ | 2 (y + > ara) | dx 


l._¢ a | Pa iy (2#s) ) . 
d + > YonA“"\ Yon = [ ay /. dx}, 4.5) 


and the first approximation to N which takes account of the convection is 


N 


o 8 
N= 7 + yA’. (4.6) 
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The coefficient 6, does not appear in the horizontal transfer of heat, but it does determine 
the first approximation to the vertical flux of heat. It is easy to verify, from Eqs. (4.3), 
that the rate of release of potential energy that is associated with the vertical heat flux 
determined by the temperature and velocity distributions 6, and y, is equal to the rate 
of viscous dissipation of energy associated with the velocity distribution y, . Vertical 
flow of heat is the driving mechanism of the convection, and horizontal flow of heat is 
an accompanying effect, which is of higher order in A. The directions of the heat fluxes 
through different parts of the vertical boundaries produced by the temperature dis- 
tributions 6, and @, are shown in Fig. 1. 

The magnitude of A for which the above power series expansion is likely to be useful 
can be estimated roughly in the following way without going through the labour of 
numerical solutions of (4.3) and (4.4). For values of //d not too different from unity, an 
approximation to the solution of V*y, = 1 is known (Love 1927, Chap. 22) to be given 
by Grashof’s formula, viz. 


9 4\-1 2 . a 
"= 2 (1 + 2) v(t - x) y(1 — y)’. (4.7) 


6, is anti-symmetrical about x = 1/2d, so that we can regard 6, as being determined by 
(4.7) and the second of equations (4.3) together with the condition 6, = 0 on the boundary 
of the rectangle y = 0, 1, x = O, 1/2d, (taking, for definiteness, the case in which the 
boundaries x = 0 and x = I/d are perfectly conducting). As expected, the above formula 
shows dy,/dx to have the same sign over this area, with a single stationary value at the 
centre. For the purpose of obtaining an estimate of the magnitude of @, , it will make 
little difference if we take the right hand side of the second of equations (4.3) as constant 
and equal to the average value of — dy,/dx, with y, given by (4.7), over the region 
0 <2 <1/2d,0 < y < 1. This average value is 


1 } - -1 
- eo R(1+5) - 
The equation for 6, has now been made identical with that describing torsion of an 
elastic rectangular prism, the solution of which is known (Love 1927, Chap. 14). The 
solution is available in series form, and the maximum value of 6, , which occurs at y = 
1/2, x = 1/Ad, is 





1 fey = 
‘ 4 4 3 ‘ 
2880 (1 + l*/d*) : == (2n + 1) cosh (2n + 1)xl 


4d 


For 1/d = 2, this maximum value is 0.97 X 107‘, and for l/d = 4 it is 0.79 K 10™. 

Thus for values of //d not very different from unity the maximum value of 6, may be 
expected to be about 10°*. When A has the standard value 1000 d*, the maximum value 
of the second term A@, in the power series for @ is then about 1/10 d*, whereas the value 
of the first term y at the same point is 1/2. It seems that even for the smallest values of 
d (and, consequently, of A) likely to be of interest (viz. d ~ 1 cm) the second term is not 
far from being as large as the first. It is fairly clear from the form of equation (4.4) that 
likewise 6, will not be a small perturbation of the temperature distribution y + AQ, at 
these same values of A. It is of course still possible that the power series (4.1) and (4.2) 
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are convergent at these same values of A, but the more important practical point is that 
the first few terms of the series—and only the first few could be determined numerically 
without excessive labour—do not by themselves provide a good approximation to the 
sums of the series. These remarks have been justified only for the case in which 1/d is 
not large compared with unity, since it is only for such values that Grashof’s formula 
(4.7) applies; however it is very unlikely that convergence of the series could vary 
appreciably with the value of l/d. It seems that A = 1000 is outside (although possibly 
only just outside) the range of usefulness of the expansion in powers of A. 

When //d > 1, Grashof’s approximate formula for y, fails because it spreads the 
variation of y, over the whole of the range 0 < x < 1/d, whereas the solution of V‘*y, = 1 
will be such that the distribution of y, with respect to y tends to a constant form as the 
distance from either end of the rectangle increases.* This asymptotic distribution of 
v, , which will surely be a valid approximation when e < x < I/d — e and eis only a few 
multiples of unity, is 


“= 4 y°(1 — y)’. (4.8) 


In the part of the rectangle where (4.8) holds, 0, , ¥2 , 6. , --+ are all zero**, correspond- 
ing to the fact that @ = y, ¥ = Ay, , is an exact solution of the full equations. Thus at 
values of A such that the series expansions (5.1) and (5.2) are rapidly convergent, finite 
values of @ — y are confined to regions at the two ends of the rectangle which have 
dimensions in the z-direction only a few times the width of the rectangle. This again 
reveals the serious limitations of the range of usefulness of the expansion in powers of A. 
Anything resembling a boundary layer on the vertical walls, or having strong asymmetry 
about y = 1/2, could not be described by the power series. In short, only at values of A 
such that conduction has a much stronger influence on the temperature distribution 
than convection will the power series be useful***, whereas measurements at practical 
value of A show that the two are of comparable importance. 

5. General value of A, andl/d-—+~. The starting point for an investigation of the 
convection at larger values of A is supplied by the remark, at the end of the last section, 
that when //d is large enough the variables 6 and y take up their asymptotic form 


A 2 a 
6=y wv=5,00-, (5.1) 


at all points not near to either end of the cavity (see Fig. 2). This asymptotic state 
corresponds to a purely vertical flow of the fluid in which the rate of viscous dissipation 
of kinetic energy, per unit length (in the vertical direction) of the cavity, is just equal 


*In accordance with Saint Venant’s principle in elasticity theory. 

** Although the difference between 6, ¥n(n > 1) and their zero asymptotic values probably increases 
as n increases, since the solutions of equations like (4.3), (4.4) ... are such that the dependent variable 
lags behind the right hand side when the latter tends to zero at infinity. 

***It is relevant that the smallest value of the Rayleigh number at which the fluid between two hori- 
zontal boundaries at different temperatures is unstable, i.e. at which convection can overcome the damp- 
ing effect of conductivity and viscosity, is also of the order of 10°; for instance if the two boundaries are 
rigid perfect conductors, the critical value of A (based on distance between the planes) is known to be 


1708. 
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to the rate at which potential energy is being released in unit length of the cavity. Since 
the horizontal velocity is zero the horizontal heat flow, over parts of the cavity where 
(5.1) holds, occurs wholly by conduction. Now at small values of A such that 

6=y+ 0,A, y = yA, 
is valid everywhere, the region of the cavity in which 6 and y do not have their asymptotic 
form (5.1) is confined, as already stated, to nearly-square regions at the two ends of the 
cavity. As A increases, and more terms in the power series become significant, the de- 


1/2 





1/2 
— 


uzSy(l-y) (I- 2y) 











Fic. 2. The asymptotic state far from both ends of the cavity. 


partures from the asymptotic state spread from each end of the cavity (especially on 
the side of flow away from the end). At large values of A at which boundary-layer 
analysis is appropriate, the asymptotic state will exist only at distances sufficiently far 
from each end for the boundary layers on the two vertical sides to have become thick 
enough to overlap and amalgamate. 

In this section we shall assume that the value of l/d is large enough for the asymptotic 
state described by (5.1) to be set up over a finite range of values of x, say X/d <x < 
(1 — X)/d, at the general value of A under consideration. Several useful conclusions 
follow from the existence of the asymptotic state at some values of 2; 

(a) The horizontal flow of heat due to convection occurs at the two ends of the 
cavity, i.e. within the ranges x < X/d,x > 1 — X)/d, and the total heat flux though a 
vertical boundary due to convection is independent of /; an increase in | merely increases 
the extent of the region in which the asymptotic state occurs and in which the transfer 
is by conduction alone. The general relation (2.11) thus reduces to 


N = 1/d + F(a, A), (5.2) 


where the function F describes the additional heat transfer due to the existence of con 
vection; likewise the addition to the conductance C due to convection is inversely 
proportional to 1. 

(b) It follows from (a) that the insertion of horizontal partitions in a deep cavity 
can only increase the total horizontal heat flow; indeed, if the new cavities which are 
created by the insertion of horizontal partitions are themselves deep enough for the 
asymptotic state to be set up in each separate cavity, we shall have 


N = 1/d + nF(o, A), 


where / is the over-all height and n the number of separate cavities. 
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(c) The streamlines which rise on the side of the cavity adjacent to the hot boundary 
ultimately become the streamlines of the falling stream on the other side of the cavity. 
In the course of its journey from one side of the cavity, in a region of asymptotic flow, 
to the other, the rising stream loses heat at a rate equal to the net flux of heat across 
any horizontal plane in the range X/d < x < (h — X)/d, i.e. equal to 


ak. ae . 
(1 — 2y) ay E y (1 — | dy 


+1 1/2 


KT, — To) | (WO) xsacre(t-xy/sa dy = kK(T, — T) [ 


A 
= KT, — To) = 5.3 
1 0) 720 ( ) 
This addition to the heat content of the air in the upper end of the cavity (i.e. in the range 
0 < x < X/d) is that due wholly to the existence of convection, and must be lost by 
conduction through the two vertical walls* in this region, some of it—a fraction a, 
say,—through the cold boundary, and the remaining fraction 1 — a@ through the hot 


boundary. Hence the contribution to N due to the existence of convection is 


A 
F(o, A) = (2a — 1)=——5 (5.4 
’ 720 ) 

a is itself a funetion of A, about which a little information is available. When A is 
very small we have already seen that the heat conveyed by the rising stream is given to 
the two vertical boundaries in nearly equal parts (because the effect of conduction is 
so powerful), the value of a being found from (4.6) and (5.4) as 


a = } + 3607.A. (5.5) 


As A increases, the distribution of loss of the heat of the rising stream becomes more 
asymmetrical and @ increases, presumably monotonically. a can never be as large as 
unity**, sothat as A © it must asymptote to a constant value, 8 say, which cannot be 
far from unity. 

We thus have the result that, provided //d is always large enough for the asymptotic 


flow to be set up for some values of x, 


l 28 — | 

Nw } is 5.6 

d 720 (5.6) 

for large A, where 8 is a constant. There is the interesting consequence that for a given 
value of 1, N (and also the conductance C) has a minimum at a value of d given by 


l (78 _ 1) 
= ane P. (5.7 
d 240 ) 
provided that this value of l/d is large enough for a region of asymptotic flow to be set 
up). With the standard conditions, and the approximation 6 = 1, this gives d = 
(1/4)'“* as the optimum spacing between the cavity boundaries, where / and d are 


*Assuming that the flow of heat through the horizontal boundary is negligible, either because its 
length is small or its conductivity is small. 

**For as the rising streamlines approach the upper end of the cavity, those near the boundary neces- 
sarily diverge and reduce the value of (06/dy),-, , so that it is less than unity, and the flux of heat through 
this part of the hot boundary into the cavity is less than that which occurs in absence of any convection. 
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expressed in centimetres. At this optimum value of d, convection is responsible for 25% 
of the total transfer, the proportion increasing very rapidly with increase of d. 

Before the significance of the results of this section can be assessed we must estimate 
the smallest value of //d for which the asymptotic state described by (5.1) is set up. It 
is difficult to do this exactly, but a number of consistent estimates can be obtained by 
different methods. There is first of all the simple argument that when the rising stream 
turns around at the top of the cavity, a new temperature distribution across this stream 
has to be set up by a process of conduction across the flow and the time needed for this 
process is of order d’/4x. During this time a particle moving with a velocity equal to the 
mean downward velocity in the asymptotic region (i.e. equal to the average value of 
Ax/12d y (1 — y) (1 — 2y) over the range 0 < y < 1/2, which is found to be A«/192d) 
moves through a distance 


A d, (5.8) 


which will be of the same order as X. 

Another estimate of X/d can be obtained by assuming that the rate at which the 
influence of the cold boundary spreads into the stream, after it has turned around at the 
upper end of the cavity, is the same as the rate at which the thermal boundary layer 
grows on an isolated cold plate by free convection. On defining the edge of the thermal 
boundary layer to be where the temperature (relative to that of the ambient fluid at 
infinity) is 20% of the temperature of the boundary—assumed to be T, — T,—we find 
from Schmidt and Beckman’s experiments (Goldstein, Ed., 1938) that the thermal 
boundary layer due to free convection from a vertical plate has a thickness d/2 when the 


length X of the plate is given by 





d stad ud _— mo) 

aa 2x 4yT, ‘ 
X A - 
ie. = = opa0' (5.9) 


Finally, we may analyse the way in which the thermal influence of the cold boundary 
spreads into the stream by means of the Graetz-Nusselt procedure (Goldstein, Ed., 
1938), which was developed for the case of forced convection due to Poiseuille flow in a 
circular tube with a sudden change in temperature of the wall at a certain section. If 
the falling stream on the cold side of the cavity is regarded as a forced flow, with parabolic 
velocity distribution (more accurate representation of the real velocity distribution is 
not warranted) and the same mass flow (for 0 < y < 1/2 only) as that given by (5.1), 
the approximate equation satisfied by @ in the falling stream is 

A a6 "6 . 

g 44 -Da = af (5.10) 
A particular solution is 

@=y+t+e™”* fly), 
provided f(y) is such that 


f” + ry — wf = 0. 
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If f(y) is now assumed to be an even power series in (y — 1/4) (the odd powers represent 
a temperature distribution which dies out, as x increases, more rapidly than that described 
by the even powers), the condition that f(y) = Oat y = 0 and y = 1/2 gives an equation 
for \ whose degree is higher as higher powers are included in the series for f. In this way 
the first three successive approximations to the smallest root for \ are found to be 615, 
735 and 705, and the series corresponding to this last approximation to the root is 


f(y) = Fil — 22.1(y — 2)? + 139y — )* — 840y — 9° + --°]. 


The general solution is then 


@=y+ Die” f,(y), (5.11) 


where \; , A, , -:: are the possible values of \, and the constants F,, are chosen to give a 
correct representation of the temperature distribution at c = 0. Whatever the shape of 
the temperature disturbance, the rate at which its magnitude diminishes as z increases 
will ultimately be dominated by the term in the series (5.11) corresponding to the smallest 
value of \. The distance over which the magnitude of this term falls to a tenth of its 


initial value is 
A 
9450 ’ (5.12) 
which can be regarded as another estimate of X/d. 

These three estimates of X/d, (5.8), (5.9) and (5.12), are consistent (some variation 
is to be expected, since X is not a definite length), despite the fact that two of them use 
the data of forced convection and one uses free convection (although the methods are 
not essentially different, since the assumed velocity of the forced convection is obtained 
from a consideration of the velocity that would ultimately be attained in free convection). 
As a rough average of the estimates, the value of z can be taken to be 

-_ © 

d 1000 ’ 
this gives the distance from either end of the cavity at which the asymptotic state 
described by (5.1) is attained approximately*. Hence for the asymptotic state to be set 


up in a cavity we must have 





l A . 

1 > 500° (5.13) 
With the standard conditions, and / and d in centimetres, this is equivalent to 

i> of, (5.14) 


showing that results obtained in this section will be relevant to many cases of cavities 
used in buildings, especially those in double windows in view of the common choice of 
the value d = 2 cm. 

*Note that it is only when A is less than about 1000 that conduction can carry the thermal effect of 
the boundaries across the flow in such a short time that X is of the same order as d; this is consistent with 


the conclusion, obtained in section 4, that A = 1000 is near the end of the range of values of A at which 
conduction is sufficiently dominant for the expansions in powers of A to be useful. 
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tw 
bo 


6. General value of //d, and A —-~. Since it will also be true, according to the 


criterion (5.14), that the value of l/d in some practical cases is not large enough for the 
asymptotic flow to be set up, it will be useful to supplement the foregoing theory with a 
consideration of the extreme case in which the rising and falling streams at the vertical 
boundaries are far from being contiguous. The thickness of boundary layers produced 
by free convection always decreases as the Rayleigh number increases, so that this 
ed by holding //d constant and allowing A to approach infinity. 


extreme case is achiey 
, the thickness of the region of thermal influence of each portion of 


Ultimately, as A 
the boundary becomes smal] compared with both / and d, and a single continuous bound- 
ary layer surrounds the cavity. The boundary layer thickness remains finite, in spite of 
the fact that the path length of any fluid particle in the boundary layer increases in- 
definitely, because gravity produces a force on the particle in the direction of flow (for 
part of the time) and counteracts the loss of momentum due to viscous forces. In the 


core of the cavity, the flow is presumably such that the influence of the diffusivities v 


and «x is slight 

This is a type of flow problem which has considerable intrinsic interest, as Pillow 
(1952) has pointed out in a discussion of the analogous situation in which heat is being 
transferred, by means of a two-dimensional convection call, between two horizontal 
plane boundaries. The special interest—and the difficulties—are associated mostly 
with the determination of conditions in the core of the cavity. Inasmuch as when A 


the terms containing second-order differentials in the governing 


e) 


(i.e. vy > 0, kx — O 
equations (2.7) and (2.8) become negligible everywhere, except possibly in the neighbour- 
hood of the boundaries, the flow in the core of the cavity conforms to 


0(6, w) 1 A(w, w) 00 
= = 0, = A : * (6.1) 
Ol7r, y) o OUN, y) OY 


Now in cases in which the flow outside a boundary layer is known to be free from vorticity 
all streamlines originate at infinity where the vorticity 


(usually by reason of the fact that 
L elocity 


is zero), this same flow is determined uniquely (by the condition that the norma 
at rigid boundaries is zero). In the above case of cavity flow we have no reason to expect 
irrotational flow in the core, and a unique determination from the condition of zero 
normal velocity at the boundary is not to be expected without the introduction of 


information concerning the flow outside the core of the cavity. Pillow shows that the 
solution of equations (6.1) is 

6 = f(y), w cAxf'(y) + gly), (6.2) 
where f and g are arbitrary functions. The condition of symmetry of the equi-tempera- 
ture lines and streamlines in the core about Z l 2d and about y= l 2 necessitate that 
f() is constant, and by symmetry again the value of the constant is 1/2. This leaves the 
vorticity distribution arbitrary, apart from being symmetrical, and to determine the 


function g(W) it would seem to be necessary to make use of the condition of compatibility 
of the flows in the core and in the boundary layer. So far as the form of gi) is concerned, 
we can employ the argument that if the vorticity in the cavity core were not uniform, 
the existence of a finite, although small, viscosity would make it uniform in a sufficiently 
long time (unless of course there is some mechanism generating vorticity in the core; 
inasmuch as the temperature is uniform in the core this possibility can be rejected). 
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Hence the conclusion is that 
6 = }, w® = w,(const.), (6.3) 


in the core of the cavity. 
The flow in the cavity thus consists of an isothermal core in which the velocity dis- 


tribution is given by 
V'v = —"@ ’ (6.4) 


with ¥ constant on a rectangular boundary, surrounded by a continuous boundary layer 
in which the temperature and velocity make rapid transitions to their assigned values 
at the boundary. If s denotes (non-dimensional) distance along the cavity boundary 
from the point A towards B in Fig. 1, the velocity outside the boundary layer can be 
written as w, U(s), where the function U is periodic in s with period 2(1/d + 1) and is 
determined by (6.4). (For the long narrow cavities with which we are concerned, U(s) 
will have the value 1/2 at points not near the two ends of the cavity, since the velocity 
will vary linearly across the cavity near such points). Except in the immediate neighbour- 
hood of the corners, the boundary layer flow will be the same as that produced on a plane 
wall, with the (spatially) periodic velocity U(s) at the wall in the absence of the boundary 
layer, periodic temperature conditions at the wall, and periodic variation of the buoyancy 
force. If u and v now denote velocity components parallel and at right angles to the 
boundary, whatever its direction, the (non-dimensional) equations describing the flow 
in the boundary layer will be, after the usual approximations are made, 


Ou Ou on, OU : au “ae 
Ua. +1 —- wo l zm oGA(@— 4) +6 ane? (6.5) 
00 00 0°60 ‘ 
Y 9s vr dn an?’ (6.6) 


where n denotes the normal distance from the boundary and G is a periodic factor which 
takes account of the change in relative direction of the buoyancy force and has the 
values shown in Fig. 3. The conditions at the wall are as before, while when n becomes 
large the flow in the cavity must be recovered, ie. 


asn >, 6—3 and u — w)U(s). (6.7 
n As now, {5~n 
ts 





D38oA = @x0 B280C os: O A 
G =0 Gs! G=0 G=-! 
Fic. 3. (The points ABCD correspond to the points ABCD in Fig. 1). 


The value of w, is presumably to be obtained mathematically from the condition 
that equations (6.5) and (6.6), with the boundary conditions specified, do in fact permit 
a periodic solution, although I have not been able to establish this in detail. Physically, 
one may say that a finite positive value of w, is produced by the torque exerted on the 
cavity core by the boundary layer surrounding it. The velocity in the boundary layer is 
everywhere in the direction of s increasing (which is obviously true immediately after 
fluid at temperature 6 = 1/2 throughout the cavity is released from rest, while in the 
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steady state a negative value of u(s) would imply the existence of a stagnation point in 
the flow which would be contrary to the whole notion of the existence of a continuous 
boundary layer) and so a core at rest would be subject to an anticlockwise frictional 
couple at all points on its perimeter. Probably the value of w. is such as to make the 
velocity at the edge of the core comparable with the maximum velocity in the boundary 
layer produced by free convection on a plate of length / and temperature difference 


1(T, — T>); that is, using the measurements of Schmidt and Beckman (Goldstein, 
Ed., 1938), 
1 wok — A faAyi/2 g(T, — my” 
2 d —_ 0.541 ! oT, , 
giving 
l 1/2 
Wo = 0.76( 0 1) : (6.8) 


I have not found any reasonably simple method of solving the boundary layer equa- 
tions (6.5) and (6.6). Methods which rely on a polynomial representation of @ and u are 
frustrated by the fact that @ and u make several oscillations as n increases and poly- 
nomials of small degree in n are quite inadequate to represent them. Linearization of 
the equations in the Oseen manner was also tried, but the simplification is evidently 
too drastic, since a solution periodic in s did not appear to be consistent with any value 
of wo. 
Even though the problem remains unsolved in detail, it is possible to estimate the 
rate of heat transfer. The flow near each vertical wall bears a partial resemblance to 
forced convection due to a stream of speed wox/(2d) past a plate of length / with tempera- 
ture difference T, — 7, (for the boundary layer near each end of the cavity, and con- 
sequently near the leading edge of each vertical wall, can be assumed to have the tempera- 
ture of the vertical wall it has just left*), on the one hand, and to free convection past a 
plate of length / with temperature difference 7, — T, on the other. In the former case 
the heat transfer would be (Goldstein, 1938), for air, and laminar flow, 


. ico i xl\'”? 
N = 0.50(4 Wo a!) s 


A 5) 
= D —_-as ).¢ 
0.38( Fu (6.9) 
in view of (6.8), whereas in the latter case the heat transfer would be (Goldstein, Ed., 
1938) 
: A ‘~ 1/4 
N= 0.48(4 3) , (6.10) 
od 


The two expressions have the same functional form (necessarily, in view of their es- 
sentially similar bases), and the closeness of the numerical coefficients makes it un- 





*It seems to be appropriate to use (7; — Jo) as the temperature difference of the equivalent isolated 
plate in an estimate of the maximum velocity in the boundary layer, as has already been done, since the 
maximum velocity occurs near the trailing edge of the plate and the temperature of the cavity core will 
here be important; but to use (7; — 7») as the temperature difference in an estimate of the heat transfer, 
because this is greatest near the leading edge where the ambient fluid is the oncoming boundary layer 


from the other vertical boundary. 
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necessary to consider which of the two cases represents the real situation more closely. 
N is now independent of d, provided //d is sufficiently large for the boundary layer to 
be dominated by its travel over the vertical walls. With the standard conditions, and 
with / in centimetres, the mean of (6.9) and (6.10) reduces to 


N = 2.01. (6.11) 


These expressions for N are asymptotically valid, as A +, and supplement the ex- 
pression (5.6) which is valid when l/d is large enough for the boundary layers to have 
amalgamated completely near the centre of the cavity at least. 

7. The criterion for the flow in the cavity to be laminar. The discussion of the 
conditions under which the flow may be expected to be either laminar or turbulent has 
been delayed until now, since it is necessary first to know what form the laminar flow 
regime would take before considering its stability. Since the parameters A and l/d 
(leaving aside o) are sufficient to determine the flow uniquely, the criterion for the flow 
to be on the borderline of the laminar state can be expressed as a relation between A and 
l/d. The kinds of laminar flow which exist when either //d or A is very large have been 
described (in sections 5 and 6), and the criteria for these two flows to break down will be 
considered sep rately. 

Taking first the case in which A is very large and a continuous boundary layer 
surrounds the cavity, we can make use again of the general resemblance which this 
boundary layer bears to the boundary layer produced by free convection on an isolated 
flat plate of length / and temperature difference 7, — 7). Experiments with air (Goldstein, 
1938) suggests that this latter flow is laminar provided the Rayleigh number based on 
the length of the plate is less than about 10°. We may expect, therefore, that an approxi- 
mate criterion for the flow (of air, and perhaps of other fluids also) in the cavity to be 
laminar at large values of A is 


3 


A LU < 10°, (7.1) 
d 


and the corresponding relation between A and //d expressing the borderline state is 
shown in Fig. 4. With the standard conditions, (7.1) states that the flow in the cavity 
will be laminar proyided 1 < 100 em (and note that the initial assumption that A is 
large here means that d must be sufficiently large). 

In the other extreme case in which l/d is so large that the asymptotic distribution of 
temperature and velocity is set up at positions not near the ends of the cavity, it seems 
likely that breakdown of the laminar flow will occur first in the region far from the cavity 
ends* since this is where the greatest velocity gradients exist. We thus require to know 
the critical Reynolds number of the flow described in Fig. 2. The boundary layer pro- 
duced by free convection on an isolated plate has a velocity distribution with similar 
general characteristics and again we may use it as a guide. The critical Reynolds number, 


formed from the maximum velocity and the boundary layer thickness, in this latter ease 


*There is also the possibility of breakdown of the laminar flow near the ends of the cavity where the 
streamlines are curved and the circulation decreases with increase of distance from the centre of curva- 
ture. A comparison with the case of flow between concentric cylinders, with the outer cylinder stationary, 
suggests that the cavity flow is near the limit of stability for some of the relevant values of d and I, but 
since there is so little time for a disturbance to be amplified while turning the corner this region of the 
flow is unlikely to become a source of turbulence. 
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is about 300 (Goldstein, 1938). Hence an approximate critical relation between A and 
1/d for the cavity flow when I/d is large, obtained by taking d as equivalent to the bound- 


ary layer thickness and by taking the maximum velocity difference in the asymptotic 
cavity flow as equivalent to the maximum velocity in the boundary layer, is 


e.-7ee Ps » = 300 
72 V3d 4% 


A = 18700c, = 13700 for air. (7.2) 


i.e. 
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Fig. 4. The critica] relation governing the existence of laminar flow. 
With the standard conditions, (7.2) states that when l/d is sufficiently large, d = 
(13.7)'”* = 2.4 em is about the largest value of d for which the flow is laminar. 

The two critical relations (7.1) and (7.2) cross at 1/d = 42**, as shown in Fig. 4; 
presumably in the neighbourhood of the point of intersection neither of the approxi- 
mations on which (7.1) and (7.2) are based is valid. The true critical relation between 
A and I|/d will have the solid portions of the two curves as asymptotes and will lie a 
little above both of them in order to have a smooth transition from one to the other. 
For practical purposes it will be sufficiently accurate to regard the flow in the cavity as 
laminar if the point (l/d, A) describing the flow lies below either of the curves (7.1) or 
(7.2). If the point lies only a little above both curves, the flow in the cavity may be 
turbulent, but it will probably be necessary for the point to lie well above the curves 


**This value of |/d is implicitly an estimate of the minimum value of //d, at A = 13700, for which the 
asymptotic velocity and temperature distributions are set up and as such is consistent with the estimate 
which is given by (5.13) (viz. 27); likewise the value A = 13700 is implicitly an estimate of the minimum 
value of A, at1/d = 42, for which the notion of a continuous boundary layer surrounding a core in which 
temperature and vorticity are uniform is valid. 
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before the turbulent region is large enough to dominate expressions for the heat transfer; 
for instance, the variation of N with Rayleigh number for the free convection produced 
by an isolated vertical plate does not settle down to the form appropriate to wholly 
turbulent flow until the Rayleigh number (based on length of the plate) is above about 
10°". As a consequence, it will be possible to use the results appropriate to laminar flow, 
as a reasonable approximation, even when the values of //d and A are not quite such as 
to satisfy the above criteria for laminar flow. 

8. Cases of turbulent flow. It seems probable, from the considerations of the previ- 
ous section, that in some cases of cavities used in buildings the convective motion is 
turbulent. The difficulties of the problem are already great when the flow is laminar, 
and it would be unwise to speculate about the case of turbulent flow without some 
guidance from experiments. All that will be done here is to indicate how the expressions 
for heat transfer obtained earlier, for the cases in which either //d or A are very large, 
are altered if the flow is assumed to be turbulent instead of laminar. 

Consider first the case described in section 5, in which I/d is so large that the two 
boundary layers on the vertical walls have amalgamated and do not vary with x over 
the range Xd < x < (l — X)/d, and assume now that the flow in this latter range is 
turbulent. The streamlines of the mean flow are vertical, as before, but there is now : 
fluctuating velocity in the direction perpendicular to the vertical boundaries and the 
convective flux of heat is not zero in this region. The magnitude of the velocity fluctua- 
tions is readily found by returning to Eqs. (2.3), (2.4) and (2.5) and taking the mean 
value of both sides. If u, v, w now represent (dimensional) velocity fluctuations and U 


is the mean velocity, these equations show that in the region of asymptotic mean flow 


> 


. @ @-T8 fr 2. 
(a) = » GE — 1 | ((8) — (8)ya12) dy, (8.1) 


where indicates the mean value of the enclosed quantity. Except in the immediate 
neighbourhood of the walls, the first term on the right hand side will be small, as in all 
turbulent flows. The general magnitudes of u and v will be equal, and roughly constant, 
over the central part of the turbulent where the influence of neither wall is dominant, 
and (8.1) thus shows that 

(u’), (v’), | (uv) |, gTs 7 T'o)d 


are all comparable in magnitude. 

Now the rate of heat transfer across the mean flow is pc,(7, — T >) (v6), per unit 
area, and the magnitude of (vé) will be the same as that of (v’)'”” since the range of varia- 
tion of @ is unity. The value of the Nusselt number N for a cavity of height / across 
which the heat is transferred at the above rate over the whole of either vertical boundary 


is then 


a ee (8.2) 


adie K(T, — T>) k d 


This is not the whole of the heat transfer, for some heat is convected upwards by the 
rising mean flow until it turns around at the upper end of the cavity and loses heat to 
the cold boundary. The rate at which heat is convected upwards across a horizontal 
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line in the region of asymptotic mean flow is 


pl/2 


pc,(T, — To)d | U6) dy, 


0 


and if we make the rough (but reliable) assumption that (u*)'” and U are of the same 


order of magnitude, this quantity has the magnitude 


(7 ma so ig r 1/2 
w,| oie (8.3) 


A certain fraction of this—not less than 1/2, not more than unity—is conducted through 
the cold boundary when the rising stream turns around at the top of the cavity, so that 
the complete expression for N is 

(8.4) 


N =), U (@A)'”? + r.(oA)'””, 


d 
where \, and X, are constants of order unity. This expression shows what might have been 
expected, that in view of the ability of turbulent fluctuations to transfer heat across the 
cavity even when the mean flow is vertical, the upward convection of heat by the mean 
flow is comparatively unimportant. (It is to be expected that A, will be a little larger 
than i, from the nature of the above approximations, but //d is large compared with 
unity and this will make the first term in (8.4) dominant.) 

In the other limiting case, ‘considered in section 6, in which A is so large that the 
flow consists of an isothermal core of uniform vorticity, surrounded by a continuous 
boundary layer, it may likewise be possible to make modifications to allow for the 
existence of turbulence. If it can be shown from experiments that the flow in the iso- 
thermal core remains laminar, still with uniform vorticity, while that in the boundary 
layer is turbulent, it will be possible to make use of the general resemblance between 
this boundary layer and that produced by free convection on an isolated plate of length 
l and temperature difference 7, — 7, or 1/2 (7, — T>) (according to the quantity under 
discussion). The empirical expression for the heat transfer, due to turbulent free con- 
vection, from an isolated plate with temperature difference 7, T, in air is (Goldstein, 
1938) 


N= 0.134" ! ; (8.5) 
d 


which may be identified with the heat transfer across the cavity as in the case of laminar 
flow. 

9. Summary of the results and comparison with experiment. A rather confusing 
picture of different types of flow occuring at different values of A and //d has been built 
up in the preceding sections, and it will be useful now to recall the principal results and 
to see how they fit together. A comparison with the available measurements (leaving 
aside those for which all the relevant data, such as the value of l/d, is not specified) 
will be made, although these are so few in number as to leave the comparison indecisive. 
The measurements to be used are those described in the book by Jakob (1949), 
of them having been made by Mull and Reiher (1930). 

The final aim of the analysis is to determine the rate of heat transfer, and the sig- 


most 
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nificance of the various types of flow can be considered best with the aid of a diagram 
showing the variation of N over a wide range of values of A. We shall use Fig. 5 for this 











Fig. 5 
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denotes multiple point 
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purpose although the representation is designed to show clearly only the variation of 
N over the range of values of A likely to be of practical interest in problems of thermal 


insulation of buildings. The ordinate, N — lI/d, is the addition to the Nusselt number 
(see (2.10)) due to the existence of convection, giving the curves for different values of 
l/d a common value at A = 0. (The contribution due to existence of convection is not 


dominant at the values of A used in the figure, so that mental allowance for the con- 
duction term 1/d must be made.) The abscissa has been chosen as A’’*, corresponding to 
a linear variation of the important parameter d when 7, — T, is kept constant; with the 
standard conditions described in section 3, and with d in centimetres, we have 


A™ = 10d. (9.1) 
Considering first values of A less than about 1000, the appropriate prediction is that 
N — l/d varies as A®* (see (4.6)), with a constant of proportionality which very roughly 


is of the order of 10“. We saw that at these low values of A at which a power series ex- 
pansion of @ and y is useful convection is much less important than conduction. Indeed 
we shall see that the values of N — l/d at even larger values of A, at which another type 
of flow occurs, are too small for their exact magnitude to be of much interest, so that 
gor pract ical purposes we can pass over the predictions about N described in section 4. 








230 G. K. BATCHELOR [Vol. XII, No. 3 
The first curve shown in Fig. 5, at low values of A, therefore, is not the quadratic 
law (4.6) but the linear law (5.6), which was derived on the assumption that //d is large 
enough for a region of flow parallel to the walls to be set up near the centre of the cavity. 
A rough criterion for this to be a valid assumption was found to be 1/d > A/500 (see 
(5.13)), so that the linear law (5.6) will be valid up to values of A which will be increasing 
with J/d. The linear law (5.6) is not definite without some assumption about the 
value of 8, which was seen to lie between 1/2 and 1 and to asymptote to a value near 
unity as 4 —o. As a tentative approximation we may assume 8 = 1, and the first 
curve shown in Fig. 5 is 
a re. . (9.2) 
d 720 
When A is about 1000, the contribution to N due to convection will be small for all 
values of /‘d large compared with unity, as remarked above. 

When A is in the neighbourhood of (30)°, = 2.7 X 10*, corresponding to d 3 em 
with the standard conditions, the requirement for the law (9.2) to be valid is that l/d 
> 54 (equivalent to 1 > 162 em with the standard conditions again), which is already 
very restrictive. At values of A too large for the criterion (5.13) to be satisfied, the flow 
does not attain the asymptotic form near the centre of the cavity and there is a tendency 
for the gradients of the temperature and velocity to be largest in the neighbourhood of 
the walls. At values of A large compared with that for which the criterion (5.13) just 
fails, the assumption of a continuous boundary layer surrounding a core of uniform 
temperature and vorticity is valid and we can use the predictions of Section 6. Two 
alternative expressions for N, (6.9) and (6.10), were obtained (the methods being such 
as to ensure the correctness of the functional form of NV although not necessarily of the 


multiplicative constant), and the mean of these expressions, viz. 


N = 0.43(4 ) = 048A () (9.3) 


Co d° 
is plotted in Fig. 5 for various values of //d. There is presumably a curve of transition 
from the law (9.2) to the appropriate member of the family (9.3), and it will be noted 
that the criterion A/500 < I/d, giving the upper limit of the range of validity of (9.2), 
is consistent with the position of the curve to which the transition must be made. Short 
dotted tangents are shown on the curve (9.2) in Fig. 5 at the place where the curve 
ceases to be valid, according to the criterion (5.13), for the value of //d shown in brackets. 
When A reaches a certain value, given by the criteria (7.1) and (7.2), steady laminar 
10, 20, 40, the 


flow ceases to be possible. For the values of 1/d employed in Fig. 5, viz. 
criterion (7.1) is appropriate, and the corresponding limiting value of A’ is shown on 
the curves in Fig. 5 as a short cross stroke. (Note that the upper limit of A for laminar 
flow does not decrease indefinitely as 1/d increases; when |/d = 42 the criterion (7.2) 
takes over, and laminar flow is possible for A'’* < 24 however large //d may be). When 
l/d = 40, the largest value of A’’® for which laminar flow is possible is smaller than 

> at which transition from the law (9.2) to the law (9.3) occurs, so that 


the values of A’ 

the curve corresponding to (9.3) with 1/d = 40 in Fig. 5 has significance only as an ap- 
proximation to the turbulent flow that will occur in practice, as explained at the end of 
Section 7. It would be useful to know where the curves, that describe the variation of 
N in turbulent flow for these larger values of 1/d, occur in Fig. 5, but speculation about 


this should perhaps await some guidance from experiments. 
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Turning now to the experimental data about N, all the measurements described in 
Jakob’s book (1949), except three referring to small values of //d and very large values 
of A, are reproduced* in Fig. 5. Where several of the measurements were nearly identical, 
they have been amalgamated into a single point, the number of separate measurements 
being indicated in brackets. At small values of A, less than about 10*, there are too few 
measurements to permit any conclusions about the validity of (9.2) and all that can be 
said is that the theory and the experiments are not inconsistent. (These smaller values 
of A are not unimportant practically, so that there is a real need for further experiments 
in this range). At larger values of A the experimental points seem to be consistent with 
a set of curves like (9.3) in form but having a numerical coefficient smaller by a factor 
of about 0.6. In other words, if (9.3) were replaced by 


3/4 
N = 03a) ; (9.4) 


the two laws (9.2) and (9.4), with appropriate transition curves, would fit the above 
data adequately. That a change in the numerical coefficient of (9.3) should be necessary 
for agreement with observation is not impossible in view of the numerical uncertainties 
involved in its derivation. Moreover it will be recalled that (9.2) was based on the notion 
of independent boundary layers on the two vertical boundaries, which will be a valid 
picture only when the boundary layer thicknesses are small compared with d. Now if 
l/d and T, — T> are kept constant, the boundary layer thickness, as a fraction of d, is 
proportional to d *’*, which is a fairly slow rate of decrease, so far as the range of values 
covered in Fig. 5 is concerned; in fact it is readily found that at A’”* = 100, the boundary 
layer on each vertical boundary near the centre has a thickness of about 0.15 d for 1/d = 
10, and varies as (I/d)'’*, which is scarcely small enough for the asymptotic picture of 
completely separate boundary layers to be valid. It is possible, therefore, that the ex- 
perimental points are approaching the curves (9.3) (or a set of curves with a numerical 
factor a little different from 0.48), and that values of A below 10° are still in the transition 
range. 

The above interpretation of the measurements is quite different from that proposed 
by Jakob (1949). Jakob plots log (Nd/l) against log A, and represents the measurements 
at values of A above 10° by means of curves of the form 


1/3 1/9 
: N = 0.065( 4) (4) (9.5) 


on the supposition that the flow is then turbulent (compare (8.6)) and the measurements 
in the range 10* < A < 10° by curves of the form 


d A 1/4 d 1/9 
-~N= 0.18(4) =F (9.6) 
l o l 
The measurements at values of A less than 10* are represented by a curve common to 
all values of //d. This is an empirical representation, and will serve as well as any other 
means of describing the data. However it should be observed that there is little evidence 
for Jakob’s assumption of turbulent flow at values of A above 10°. The criterion for 
*I have been obliged to take the data from Fig. 25-7 of Jakob’s book since Mull and Reiher’s original 


paper is not available in Cambridge. 
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laminar flow to be possible clearly depends on |/d as well as on A, and the considerations 
of section 7 suggest that all the measurements reproduced in Fig. 5 (and also the three 
measurements not shown) correspond to laminar flow. The empirical conclusion that 
log dN /l is independent of //d at values of A below about 10° is also at variance with the 
theory given herein, since this is the range in which we expect a region of parallel, vertical, 
flow to be set up near the centre of the cavity, with (9.2) as the corresponding expression 
for V. 

For many practical purposes, the description of the theoretical results that is given 
in Fig. 5 is not the most revealing. The quantity of greatest interest is the heat transfer 
per unit area of a vertical boundary, which is proportional to N/l. Moreover, the width 
d can be chosen much more freely than the height / when a cavity is being used in a build- 
curves showing the heat transfer as a function of d for various values of l, 


ing, so that 
are needed. Both of these needs are met by Fig. 6, 


rather than for various values of l/d, 
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as a function of d (i.e. of 


which shows the thermal conductance C (defined by (2.12) 
A’”?. when 7, — T, is given) for various values of 1, with the standard conditions. The 


range of values of d is | 
for C are those in common use by heating engineers. 

At values of d less than about 1 em, the conductance has approximately the value 
ignoring the effect of convection, but at larger values of d the curves obtained 


ere restricted to those of greatest practical interest. The units 


obtained by 
from (9.2) with different values of | diverge as shown in the figure. (There is the useful 
conclusion that if it is possible to subdivide a given rectangular cavity by the insertion 
small thickness, a maximum width of | cm for any of the new 
cavities so formed will ensure that the overall conductance is close to the optimum 
value obtained in the absence of convection.) The curves for different values of / all 
show the potentially important feature of a minimum value of C. However the position 


of vertical partitions of 
| 
l 
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of the curves corresponding to (9.3), towards which N asymptotes as d ©, is such 
that the value of N rises only a little, if at all, as d increases above the value at which 
the minimum of (9.2) occurs. The limiting curves given by (9.3) are shown on the right 
side of Fig. 6, but no transition curves to join the two theoretical predictions have been 
drawn in view of the uncertainty about the rate at which the limiting curves are ap- 
proached. 
‘asurements described by Jakob have been plotted in Fig. 6 by the device of 
calculating what the values of d and / would have been if g(7, — TT >)/To«v had the 
standard value 1000 cm “ and if A and //d had the values given by Jakob; these equiva- 
lent values of | are shown in brackets by each point. As already seen, they are consistent 
with the theory only if we change the factor 0.48 in (9.3) to 0.3, or if we suppose that the 
law (9.3) becomes asymptotically valid at values of d well above those used in Fig. 6. 
Neither the theory no the measurements support the existence of a definite minimum 
of Cas a function of d but both suggest that no further significant decrease in C occurs 
for values of d above about 2.5 em (not even when | is greater than 200 cm, because for 
such cavities laminar flow becomes impossible when d > 2.4 cm, and this will prevent 


The me 


any decrease in (). 
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CORRECTIONS IN HOT-WIRE CORRELATION MEASUREMENTS* 


BY 
CARL E. PEARSON 


Harvard University 


Summary. Particularly in small-scale turbulence, hot-wire correlation measure- 
ments must be corrected for the effect of finite wire length. End effects, imperfect velocity 
correlation along the length of the wire and transient longitudinal heat conduction must 
be considered. An integral equation for these corrections is obtained; with appropriate 
experimental procedure, this equation may be solved explicitly. 

1. Introduction. A hot-wire instrument consists of a wire filament (e.g. tungsten 
of .005 mm diameter and 1 mm length) heated by an electric current and cooled by 
the flow of the fluid in which it is immersed. The rate of heat loss of such a wire is found 
experimentally to be proportional to the temperature difference between wire and 
stream multiplied by an empirical function of stream velocity. When the stream velocity 
is fluctuating, as in turbulence, the relationship between instantaneous values of tem- 
perature, velocity, and rate of heat loss is assumed to be the same as that between their 
steady-state values; this assumption, justified by an order-of-magnitude heat transfer 
calculation and also indirectly by the results of its use, allows a study of turbulence to 
be made. 

It is customary to maintain constant either the temperature or the current in the 
wire. In the latter case, the response of the instrument to a velocity fluctuation is modi- 
fied by its own heat capacity; for convenience, a compensating circuit to correct for this 
effect is often included. Since such circuits result in a decrease in signal-noise ratio, 
they are occasionally omitted [1]. A constant temperature instrument does not require 
compensation, and moreover has an inherently faster ‘response-completion” time; 
however, the required electronic feed-back circuit decreases the signal-noise ratio. 
Although each instrument is superior to the other in specific applications, the constant- 
current type is the more widely used, largely because of its simplicity. 

In practice, a number of difficulties arise. Among these are the fact that the two 
ends of the wire are both at stream temperature because of the large amount of well- 
cooled metal in the solder joint at each end (Fig. 1); this end effect affects the response 
of the instrument—particularly if the wire is sufficiently short to respond well to small- 
scale turbulence. The magnitude of the end-effect has been calculated by Betchov [2] 
for the special case in which all elements of the wire length experience the same in- 
stantaneous velocity fluctuation (in fact, a harmonic fluctuation). Betchov has also 
discussed the effect of possible non-linearity in the dependence of heat-transfer rate 
on the temperature difference between wire and fluid; however, not only is the existence 
of appreciable non-linearity doubtful (for example, certain experiments performed on 
thin tungsten wires by the present author did not show any discernable non-linearity), 
but it appears in any case that such non-linearity effect would be completely dominated 
by end-effects for such short wires as are contemplated herein. It may be remarked that 
even if non-linearity is included, the perturbation equation corresponding to Eq. (3) is 
still linear and precisely the same Fourier Transform techniques may be used; the algebra 


*Received June 12, 1953. 
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is merely somewhat more complicated (involving elliptic instead of hyperbolic functions). 

A second problem arises from the effect of imperfect velocity correlation along the 
length of the wire. Skramstad [4] has calculated the correction to be applied to measured 
correlation coefficients in the case where (1) there are no end-effects (2) the effect of 
transient heat conductivity along the length of the wire is neglected. It has been pointed 
out by various authors, e.g. Ref. [3], that Skramstad’s equation may be obtained by 
merely manipulating the order of integration of a certain multiple integral; some general 
remarks on such situations have been made by Uberoi and Kovasznay [5] and by Liep- 
mann [6]. This latter author includes some remarks on combined space-time integral 
products (Eq. III-17 of Ref. [6])*. Our present purpose is to derive an integral equation 
for the true correlation function in terms of the apparent function as measured by two 
identical wires; it will be found necessary to measure the correlation for a sequence of 
different values of wire separation, and with different values of time displacement. It 
will turn out that there is a particular way of measuring these correlations which allows 
the integral equation to be explicitly inverted. Only one component of the correlation 
tensor is discussed; if the turbulence has isotropy as well as the assumed homogeneity, 
and if the fluid is essentially incompressible, then it is well known that a knowledge of 
this one component determines all others. It is assumed in the derivation that the two 
hot wires used are not sufficiently close to one another to result in mutual interference. 
Although the weak signals encountered in small-scale turbulence work make the constant- 
current instrument with its inherently superior signal-noise ratio the natural choice in 
such work, the constant temperature case will also be discussed. 

2. Transform of energy equation. Consider a long thin wire supported at each end 
as in Fig. 1. The steady component of fluid velocity is perpendicular to the wire axis, 


=e 





ix 








and insofar as turbulence is concerned, the wire will respond only to that component of 
velocity fluctuation which is parallel to this steady velocity direction. At any instant, 
the requirement of energy conservation applied to a wire element of length dz yields: 


2 2m 
ce 4 as ax( 27) = (T — T,) dx f(V) + gAc ao(27), (1) 
A Ox ot 

where x (cm) is measured as in Fig. 1, A (cm’) is the cross-sectional area of the wire, 
I (amps) the wire current, r (ohm cm) the resistivity of the wire metal, k (joules/cm 
sec °C.) the thermal conductivity of the wire metal, 7’ and 7, (°C.) the temperature of 
the wire element and stream, respectively (the latter is assumed constant), V the fluid 
velocity, f(V) (in joules/em sec °C.) the heat transfer coefficient from wire to stream 


*The author is indebted to the editors of the Quarterly for pointing out this reference. 

















1954] HOT-WIRE CORRELATION MEASUREMENTS 237 


per unit length of wire, g (gm/em*) the density of the wire metal, c (joules/gm °C.) the 
specific heat of the wire metal, and ¢ (sec) the time. The resistivity r at temperature T' 
is related to its value r, at 7, by means of a (/°C.), the temperature coefficient of re- 
sistivity: 
r=r,(1 +a(T — T,)]. 

For drawn tungsten wire, some typical values are: c = .142, k = 1.67, a = .0037, 
v. 8 X 10°° at 20 °C, g = 19. For such a wire of A = 6.15 X 10°’, f(V) = .00254, 
.00362, and .00519 at 40, 100, and 250 ft/sec air flow respectively. 

The velocity and resistivity may each be replaced by the sum of a steady-state and 


a fluctuating component: 
V = V+o(z2, 0), r = R(x) + r(a, 2), 


where R(x) is the steady-state solution of Eq. (1), corresponding to d7'/at = 0: 


rf(V) | T’r,« cosh (I | (2) 


R(a) = “oe rior 


~ Af(V) cosh (P,L/2) 
where 
p? _ (V) — Praa/A 
; kA f 
Then a perturbation in Eq. (1) (J is held constant; the case of constant temperature 


will be discussed subsequently) yields: 








dr pep Or _ pl, _ cosh (Piz) | 
Ox” Pr — P, at Bal cosh (P,L/2) |’ (3) 
where 
P= pa —LtealOf/aV) _ 
6 k? kA*f(V) — P’riakA’® 


The boundary conditions are that r = 0 at z = + L/2. Taking the Fourier Transform 
in time of Eq. (3) gives an ordinary differential equation 

d’r’ - : cosh (P,2) | 

—3 — (P; A)r’ = B’| 1 -—- > 
me etre oo E cosh (P,L/2) |’ 


where 


a @ 


r'(z,r) = | r(x, t) exp (—7Ad) dt, v'(z, A) = | v(x, t) exp (—7Ad) dl, 
which is easily solved for r’ (boundary conditions unchanged). The transform of pertur- 
bation voltage e across the wire may then be calculated by 


L/2 
e’(\) = / Ir’ dx/A, 


-L/2 
to give, finally, 
BI p” v( cosh (Px) ) 


O) = AP + Pd J 1s osh (P,L/2) 


~ eosh (P,L/2) 
(- _eosh [(P} + iP2d)'z] 








cosh [(P? + iP.d)'?L V) dz =) 
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The value of e may be obtained by inversion as 


e(t) = | dx de dd F(x, d)v(x, €) exp [2A(t — 8], (5) 
where 
F(z, +) = 0 for zi > L/2, 
BI ( cosh (P;2) )( cosh [(P} + iP.d)'*2x] 
A) = ; - —-—-—— ] —] + 7 Var ee 
Plz, d 2r A(P; + tP2X) cosh (P,L/2) '? cosh [(P} + iP.d)'°L/2]} 


for \2|< Z/2. (6) 
3. Correlation measurements. By a change of variable of integration, Eq. (5) may 
be written 


ei) = | dx de dd F(x, d)v(a, t + ©) exp (—1Ae). (7) 


v—« 


Suppose now that two identical wires are arranged as shown in Fig. 2 with the main 
stream velocity perpendicular to their common plane. The midpoints of the two wires 


"2 





= 7 
r? 

x, 
Lu 


are displaced horizontally by a distance M and vertically by a distance P; let x, and zx, 
be distances measured in the same direction from the midpoints along the first and 
second wires respectively. Define the correlation function C (a, b) to be the time average 
of the product of the velocity perturbation (i.e., the component parallel to main stream 
velocity) at one point in space with that at a point distant a in space (in some direction 
perpendicular to the main stream velocity) and b in time from the first point. Then for 
points x, , x, on the two wires, 


Cia, b) = C{[M”* + (P+ is: = z,)"]' ie b} 
The turbulence field is assumed homogeneous. Now the measured correlation S for a 
> 


fixed time interval of N seconds is given by the time average of the product of the voltages 


across the wires, one measured WN seconds later than the other: 


S (M, P, N) = average of [e, (t) e. (t + N)]. 
Using Eq. (5), this becomes 
S(M,P,N) = | dx, datz dey dep dd dr, F(x, , ds) F (a2 , 2 


- C{[M* + (P + a2 — 2,)7]'", N + @& — «} exp (—”rye, — Aye). 
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Setting «, — «, = eand xz, — x, = x and making use of the fact that 


ro) 


[ de, exp (—7r16 _ 1Ar€>) = 2rA(r, + do), 
[where A is the delta function of the argument (A, + Az.)] gives 
S(M, P, N) = 2r [ dx, dx dd» de F(x; , —X2) F(a, + 2, do) 


- C{[(M? + (P + a)’]'”, N + €} exp (—a”r2). 


Taking the transform of each side gives 


[ dP dN S(M, P, N) exp (taP + iBN) = 2x [ dP dN dz, dx de dd. F(x, , — 2) 


- F(x, + x, )C{[M*? + (P + 2x)"]'’, N + e} exp (taP + iBN — rr.) 


ai [ dx, de Dy UP + 2) d(x, + 2) 


“-@ 


- a(N + F(x, , —d2) F(a, + x, \2)C{[M? + (P + 2)*]'", N + 6} 
- exp [ia(P + x) — ta(x, + 2) + tax, + iB(N + ©) — iBe — irz€] 


= ir [ dis dp C([M* + w’]'”’, ¢) exp (taw + ise) 


‘ 1 [ dx dy F(x, 8)F(y, —B) exp (tax — ia) 


“-—-@ 


whence 
| dw de C([M? + w]'”, ¢) exp (iaw + 186) 


/ ” aP dN S(M, P, N) exp (iaP + i8N) 





a co 


4x’ | | dx F(x, 8) exp (iazx) 
| e 


2 





Inversion yields 


C({M? + d’]'”, s) = [ dP dN da dg Ra J iE = 9 + OS - s)] (8) 


7 Sx | [ dx F(x, 8) exp (iax) | 
| #©—L/2 





Since F(x, 8) is known from Eq. (6), Eq. (8) provides an explicit expression for the 
true correlation function C in terms of the measured correlation S. No compensation 
circuits are contemplated in these measurements. It is worth emphasizing that for fixed 
M a sequence of measurements for different values of P are necessary; this is the device 
which allows an explicit solution to be written. The integration with respect to z, a, 
and 8 may be performed explicitly if desired; Eq. (8) then has the form 


ao 


C((M? + @)'?, 8) = | dP dN S(M, P, N)W(P - a, N - 8), 


“0 
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where the function W is probably best handled by tabulation of values as calculated 
from Eq. (8). 

Note also that S could be written as a function of the distance between the mid- 
points of the two parallel wires and of the angle between their axes and the line joining 
their midpoints. 

4. Constant temperature case. If the total resistance of the wire is maintained 
constant by means of a feedback circuit, then (letting r again denote the perturbation 
in resistivity 


[ ; raés/A = | (=) a - 


I 


and, if p is the perturbation in power input to the wire, 


= Rix) —R Of ' kA f Or’) 4 
(j) = ha On, ae > — (9 
P | L/2 Ra (3) ys ae Ra | (2) L/2 es | ) 


The last term may be calculated by writing the non-integrated perturbation equation 


pRR Qa a°r P? P or (‘“ —_ R,)(of 2) 10) 
; , —- "ros ia2a- 8 . v, 
R,A‘k Ox” ol kA 
where R, is the total (constant) wire resistance. A time transform yields 
RR a Oo r’(2, r) > (R = R ofr Ov) 
——— »/(x, X) 4+ a — (P; + iP.d)r’(az, X) = ‘ v(x. X 
R,A’*k ' ax” iia kA 
Insertion of p’ from the transform of Eq. (10) gives a certain ordinary differential 
equation whose explicit solution may be written down subject to the conditions that 
r’ (—L/2, X r’ (L/2, \) = 0. The equation may be differentiated and (dr’/dxr),,. — 


calculated from the result (it will appear on both sides). The expression 


(dr’/dx)_, 
for this quantity may be inverted and substituted into Eq. (10) which now has the form 
of Eq. (5); the subsequent analysis is analogous. 

There is little point in reproducing this latter analysis in detail because of the fact 
that the situation of present interest—viz., small-scale turbulence—is precisely that in 
which the high noise level of the constant temperature instrument makes it rather un- 


suitable in comparison with (uncompensated) constant current equipment. 
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VIBRATIONS OF TWISTED BEAMS* 


BY 
R. C. DiPRIMAt AND G. H. HANDELMAN 
Carnegie Institute of Technology 

1. Equations of equilibrium. In recent years the determination of the natural 
frequencies of twisted, cantilevered beams has received much attention. The analysis 
of a twisted beam possesses certain added complications over that of an untwisted beam 
due to the geometry of the structure. W. Prager**[1] has given a very compact vector 
representation of the equations of static equilibrium of curved, twisted beams. His 
results will be extended in this paper to include dynamic effects. Since Prager’s analysis 
is not readily available, a short summary of his approach to the problem will be given. 

A beam which is twisted but not curved can be described in terms of a straight line, 
called the center line, which is the locus of the centroids of the cross-sectional planes 
taken normal to the line. A cross-section of the beam is specified by means of the arc 
length s measured in a positive sense along the center line from a fixed origin 0. With 
each cross section, a right handed triad of orthogonal unit vectors i, j, k can be associated 


as indicated in Fig. 1. The vector i is the unit tangent vector of the center line at the 
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centroid P of the cross section under consideration and points in the direction of in- 
creasing s. The unit vectors j and k have the directions of the principal axes of inertia 
of the cross section. In the case of an untwisted beam the vectors i, j, k have fixed direc- 
tions. For a twisted beam the triad rotates about ias the point P moves along the center 
line. It is this change in orientation of the principal axes of the cross section as we move 
from point to point which creates the difficulty in the problem. 

Since the triad i, j, k moves as a rigid body, its instantaneous motion can be de- 
composed into the translation of P plus a rotation about P. If the rate of rotation per 
unit length is characterized by the vector t, we note that ¢ is proportional to i. Hence 


<= Ti, (1.1) 


*Received July 16, 1953. Work on this problem has been sponsored by the U. 8. Air Force. 

+The material in this paper was submitted as part of a thesis by the first author (now at Massachu- 
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of the requirements for the degree of Doctor of Philosophy. 

**Numbers in square brackets refer to the list of References given at the end of the paper. 
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L. 


where the scalar 7 is called the natural twist of the beam. It can then be shown that 
di/ds = 0. dj ds = rk. dk/ds —7j. (1.2) 


Let f(s) and c(s) denote the vector intensities at P of the distributed forces and 
couples applied to the beam. Let the stress resultant R(s), a force applied at the centroid 
P of the cross section, and M(s), a couple, be equipollent to the stresses exerted on this 
cross section by the portion of the beam lying to the side of increasing values of s. If 
no concentrated forces or couples are applied, the force and moment equilibrium equa 
tions take the form 


IR(s 
om + fs) = 0, (1.3) 
as 


dM(s) 


7 + c(s) +i X R(s) = O. (1.4) 
ds 


Under the action of the loads f and c the elements of the beam will undergo certain 
displacements. Points lying in the plane cross section at P of the undeformed beam 
will no longer lie in the same plane, but will form a curved surface called the warped 
cross section. In the theory of thin beams it is customary to exclude shearing stresses of 
the type which would produce a change in the right angle between two dine elements 
at P in the directions of j and k. Consequently, it is possible to define a new triad i*, 
j*, k* in which j* and k* point along the displaced principal axes of inertia at P*, the 
new position of P. Let u(s) denote the displacement of P. In order to make the triad at 
P coincide with the new triad at P*, it is necessary to add to u(s) a small rotation @(s) 
about the point P*. If g(s) and h(s) denote the linear and angular distortions of the 
beam, respectively, then the following relations are valid: 


ih (S) P ~ 
g(s) = = — @(s) X 1, (1.5) 
as 
10(s) 
a « (1.6) 
ds 


For elastic beams it is assumed that these distortions can be represented as the 
products of the corresponding stress resultants and certain dyadics A and B which 
depend on the shape of the cross section and on the elastic properties of the material 
of the beam. This is analogous to the usual stress-strain law in elasticity. More precisely 
we assume 

g(s) = A-R(s), (1.7) 
h(s) = B-M(s). (1.8) 


Moreover, since the j and k axes coincide with the principal axes of inertia of the beam, 
the dyadics A and B are in diagonal form. If we write 
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the entries 8, and 8, are equal to 1/EJ, and 1/EI, where E denotes Young’s modulus 
and J, and J, denote the moments of inertia of the cross section about the vectors j and 
k. The coefficient 8, is the reciprocal of the torsional rigidity of the beam. 

Throughout this work it will be assumed that the bending due to shear forces can 
be neglected. This is a usual assumption in beam theory, and has been discussed by 8S. 
Timoshenko for an untwisted beam [2]. Hence the linear distortions g(s) can be neglected; 
and consequently it will not be necessary to consider a» , a, , @ . These effects can, in 
principle, be included in the work to be described; but the resulting system of equations, 
which would include 6, and @, as dependent variables, would be considerably more 
complicated. 

2. Equations of motion. If the beam is assumed to undergo small transverse vibra- 
tions only, then the longitudinal motion can be neglected. Using the subscripts 0, 1, 2 
to denote components taken in the direction of i, j, k respectively, the previous assump- 
tion implies that wu) = 0. Since there are no torsional oscillations, @. = 0. Furthermore, 
it is customary to neglect the effect of rotational inertia; hence c = 0. Again we point 
out that this approximation can be removed if desired. 

Equations (1.3) and (1.4) can be written as 


aR 
+ t= 0, (2.1) 
ow +ixR=0, (2.2) 


where we write R for R(s), M for M(s), and similarly for the other vector quantities. 
Since we are neglecting the distortion due to the shear forces, Eq. (1.5) reduces to 


a 6 X i; (2.3) 
Os 


and, combining Eqs. (1.6) and (1.8), we have 
00 _ 


= . 9 
B-M. (2.4) 


Eliminating R between Eqs. (2.1) and (2.2) we obtain 
aM si, . 
—z -ixf=0. (2.5) 
Os 

Also differentiating Eq. (2.3) with respect to s and substituting B-M for 00/ds we have 


du 
ds” 


= (B-M) Xi, 


from which M can be written as 
au 
M = B":-4i X =3?. 2.6 
os” (2.6) 
It should be noted that B~' is now a two dimensional dyadic. For the case of pure trans- 
verse vibrations the distributed load or reversed effective force is 


d'u (2.7) 


f= —mA ar 
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where m is the mass per unit volume and A is the area of the cross section. Combining 
all of these results, we obtain the following equation for the displacement vector u(s, ¢): 


a E aul ; ou . 
~3 {B-i X =37 + mAi X =e = 0. (2.8) 
Os Os’) ot 

We shall consider a beam which is elastically constrained at the end s = 0, and is 
free at the end s = 1. The boundary conditions at the constrained end (s = 0) may be 


expressed in the following way: 


ixu=0 
(2.9) 
0 P 
—~ {ixu} —e-M=0 
Os 
where ¢ is a dyadic with positive entries, 
€; 0 
£ = 
0 €> 
The first of Eqs. (2.9) simply states that u, = u, = 0 at s = 0. The second, when written 
out in component form, yields 
Ou 
— — ¢M, = 0, 
Os 


or 
Ou, Jou OU> Or ° 
— — ¢ KI§— — 27 — —- — uw — ru? = 0, 
Os Os Os Os ‘ 
and 
Ou, 
—— + ¢ M, = 0, 
0s 
or 


= 


Os os” Os Os 


—_ — B14? et + 2 i + oy ul — ru} = Q. 
Physically, this boundary condition implies that the angle of deflection about the j axis 
is proportional to the moment about that axis, and similarly for the k axis. If «, = 0, 
then du,/ds = 0 and hence there is perfect clamping about the j axis. If «. = o, the 
beam is simply supported about the j axis. Similar statements hold concerning «&. 
Various combinations of «, and e, may be taken. 

The conditions at s = I are the usual conditions that there be no moment or force 
transmitted across the free end of the bar. These may be written as 


ats = I, or 
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If u(s, ¢) is assumed to be of the form 


u(s, t) = v(s)e’, 


Eq. (2.9) becomes 


- {Bi « at — mA)’i X v = 0, (2.10) 
with the boundary conditions 
ixv=0 : axy—-eB'-Saxy <0 (2.11) 
: ds * ds* - 
at s = 0, and 
7 or 
B'GxXy=0, 5 {B 7a i x w} = 0 (2.12) 


at s = /. It will be shown later that Eqs. (2.12) also appear as natural boundary con- 
ditions arising from the corresponding variational principle. 
In order to obtain a similarity principle, we introduce the dimensionless variables 


z= 3/l, w = v/l, 6 = rl. 


Then Eq. (2.10) becomes 


£ {Bi x a — mAl'»’i X w = 0, (2.13) 
with the boundary conditions 
ixw=0, 20x) —e-BY- Sa xw <0 (2.14) 
dx dx 
at x = 0, where e’ = ¢/l; and 
B os (i X w) = 0, a pS (i Xx w} =0 (2.15) 
at x = 1. If the dimensions of the cross section of a cantilevered beam [e’ = 0] are 


changed by a scale factor a, then the area A’, and moments of inertia J{ , 7; of the new 
beam C’ are related to the old ones by A’ = a7A, I{ = a‘], , Ii = aI, , l’ = al. If, in 
addition, E’ = bE, and m’ = cm then we have for the beam C’ 

a a oa d’w a’c 2 

:{B i xX 5? - 4S mAl‘r”i X w = 0, 

dx dx b 
and the homogeneous boundary conditions remain the same. From this we see that if 
\ is the natural frequency of the original beam C then 


L ba 
a c 


is the natural frequency of C’. This is exactly the same law of similarity that could be 


used for an untwisted cantilevered beam. 
The close analogy between the theory presented here and the usual theory of un- 
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twisted beams can be brought out more clearly if we set 
y=ixXw. 


Then Eqs. (2.13), (2.14) and (2.15) become 


s. B SF — mAl'd’y = 0, 


dx” dx” 
with 
y = OU, 
dy i , dy 
—e’-B -—5 = 0 
dx ' dx 
at x = 0, and atz = 1 


reise. 2 (B 2) = 0. 


dx” dx dx” 


It must be pointed out that, although these equations appear to be simply generaliza- 
tions of the differential equation and boundary conditions of an untwisted beam, the 
effect of the twist is actually hidden in the differentiation with respect to x. The two 
scalar equations for w, and w, corresponding to the vector Eq. (2.13) are 


d (EI.p) — 26 ¢ (EI,q) — qs pr ) — S(EI.)p — mAl'\’*w, = 0. 
dx* P 0° de NY dz * id hei = wa | 
r (2.16) 
d’ m — d = . dé . 2/7 4,2 | 
s (EI,q) + 26 (EI,p) + (EI,p) — &(EI,q) — mATX w, = 0, | 
da dx dx 
where 


d’w, ox dW: dé és 
p= 5 — 26 - We — OW, 
dx dx dx 


2 ‘ 
d*w, . du, dé ve 
( 3 + 26 - WwW, — OW.. 
t dx” dx dz * ; 
It is interesting to note that for a uniform beam of constant twist Eqs. (2.16) reduce 
to two simultaneous equations with constant coefficients. They are 


d‘ y ; i? ) ‘ . ‘ o 
EI, 1 — 28°(EI, + 2EI,) <4! + El,s'w, — 26(EI, + EI,) £2 
dx dx dx 
+ 25°(EI, + EI,) an = mAl'd*u, , 
EI, Hs — 25°(EI, + 2EI,) ows + EI,5*w, + 26(EI, + EI.) ¢ = 
dx dz dx 
— 26°(EI, + Ely)! = mAtNuy . 


In principle these equations can be solved exactly. 
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3. Variational principles. It will be shown that the eigenvalues \* and the corre- 
sponding eigenvectors v of Eqs. (2.10), (2.11) and (2.12) may be obtained from a minimum 
principle. Essentially we are looking for a statement of the Rayleigh type which relates 
the square of the natural frequency to the ratio of the potential and kinetic energy. 
The potential energy will consist of two parts; the first term is the usual energy due 
to bending and to this must be added the contribution of the elastic constraint at the 
clamped end. A little investigation will show that the following vector functional plays 


the role of the potential energy 


hse Behn Me - 04 “ts “*) 
D(v) = | ( x "Y).p (i < zs) ds + | (s coe ds) '* pix ds) |,’ (3.1) 


and the functional 


vl 
H(v) = | mA(i X v)-(i X v) ds, (3.2) 
the role of the kinetic energy. The subscript 0 denotes the value of the quantity in 
brackets at s = 0. 
In addition, we introduce the mixed forms 


an & d’v 7 d’u a ag d’v ie os du 


and 


al 


H(v, u) = | mA(i X v)-(i X u) ds. (3.2a) 


Since B™‘ and e are symmetric, positive dyadics it follows immediately that D(v, u) 
and H(v, u) are symmetric functionals. Also D(v) and H(v) are positive for real vectors 
v. It can be readily seen that 


Div + u) = Div) + 2D(v, u) + D(u), (3.1b) 
and 
Hiv + u) = H(v) + 2H(v, u) + Ay). (3.2b) 


The vector v(s) will be called an admissible vector if it is continuous and has con- 
tinuous derivatives of at least the fourth order, and also satisfies the constraining boun- 
dary conditions Eqs. (2.11). The lowest eigenvalue of Eq. (2.10) subject to the boundary 
conditions Eqs. (2.12) can then be found from the following minimum principle. 

Of all the admissible vectors v(s), that one which minimizes D(v) under the side 
condition H(v) = 1 is an eigenvector of the differential equation (2.10) subject to the 
boundary conditions (2.12). The minimum value of D(v) is the corresponding eigenvalue. 

It should be noted that this principle not only characterizes the lowest eigenvalue, 
but also states that the minimizing vector will satisfy automatically the free end con- 
ditions (2.12) as natural boundary conditions. The proof of this minimum principle 
will follow as a special case of the general minimum principle for the nth eigenvalue, 
which will be stated and proved later in this section. 

An orthogonality relationship between the eigenvectors of Eq. (2.10) from which 
the reality of the eigenvalues follows, can be established by means of a standard technique. 
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Let v(s) and u(s) be two eigenfunctions of Eq. (2.10), and \’ and yw’, where \* ¥ yu’, the 
respective eigenvalues. Then 


( d’v\ 


d’ l . 9 
72/8 "#3 & 752 | — mANi Xv = 0, 
(iS . ads J 


a be ae du\ 


ds* | ds*| ~ mApi Xu = 0. 


Multiplying the first equation by i X u, the second by i X v, subtracting and inte- 


grating from 0 to / yields 


- d’ , d d° - 
| { iE: ix ?¥].G x - E ix lax wha 


Jo \ds” L ds a ds” as 


—(X — p) | mA(i X v)-(i X wu) ds = 0. 
The first integral is seen to be zero, if it is integrated twice by parts; thus 


(° — pw )H(v, u) = 0. 


Since x" = ue the follow ing orthogonality relation holds: 


H(v,u) = 0. (3.5) 
This may also be written as 
Div, u) = 0, (3.4) 


which follows from the fact that 


ds” 


0 


q . d° dé 
| | (B ix Y) |x) as = r? | mA(i X v):(i X wu) ds. 


The right hand side is zero, and integrating the left hand side twice by parts and using 
the boundary conditions, Eqs. (2.11) and (2.12) gives 


| (ix <4) -B (i x T) ds +  (s ix TY)..-(B ix | = 0 
which is the desired result. : 

A computation completely parallel to that just given also shows that if v is an eigen- 
vector of the differential equation and d* is the corresponding eigenvalue, then \” = 
D(v)/H(v). Combining this result with the minimum principle, we see that the minimum 
principle actually characterizes the lowest eigenvalue. =~ 

Now assume that ” is complex; that is, \7 = o + ip, [i = W—1], where o and p 
are real, and p + 0. Corresponding to \° there will be a complex eigenvector v = r + it 
where r and t are real vector functions and t # 0. Since Eqs. (2.10), (2.11) and (2.12) 


are linear with real coefficients, it follows that A* = o — 7p is also an eigenvalue with 
corresponding eigenvector vy = r — it. If \* ¥ »’, then the orthogonality condition 
implies 

al 


| mA(i X v)-(i X v) ds = 0. 


0 














“rs 
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Now 


| | mA(i X v):(i X v) ds = | mA(—v.j + 0,k)-(—v.j + v,k) ds 


« O 
al 
= | mA{|v, | + |v. |"} ds. 
“0 
This integral cannot vanish, consequently 
vw —-rv=0 (3.5) 
which implies 
Zip = 0 
or that p = 0. Thus the eigenvalues are real and positive. 
The general variational principle may be stated as follows. Of all the admissible 
vectors v, that one which minimizes D(v) under the side condition 


H(v) = 1 
and 
Hy, a) = @, += 1,2,---,n-—1, 


where u’ is the ith eigenvector of Eq. (2.10), is an eigenvector of the differential equation 
(2.10) subject to the boundary conditions (2.12). The minimum value of D(v) under 
these side conditions on v is then the corresponding eigenvalue \% and AZ = AZ_, = 

A 

In order to prove this statement, let v be the minimizing vector and \; the correspond- 
ing minimum value. Then 

Div) = NH(v). 

Let r be an arbitrary admissible vector that satisfies the condition H(u', r) = 0 for 
i 1,2,---,n — 1. Let a be an arbitrary parameter. Then since v minimizes D(v) /H(v) 


it re that 
D(v + ar) = \2H(v + ar). 
This can be written as 
i Div N2H(v)} + 2a{ Div, r) — AZH(v, r)} + a®{ Dit) — NH} = O. 
rhe fact that this expression assumes a minimum value for a = 0 implies that 
D(v,r) — NH(v, 4) = 0. (3.6) 


Now r may be chosen as 
n—1 
r=t+ Doan’, 
i=l 
where t is an arbitrary admissible vector and the a; are determined by the condition that 


f(r, wu’) = Hit, u’) + a;H(u’) = 0, 
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which yields 
a; = —H(t, w’). 


Substituting for r in Eq. (3.6) we have 
n-1 
D(v, t) — 2H(v, t) + >> a,[D(v, u') — NA(v, u')] = 0. 
1=1 


By hypothesis H(v, u' 0; and following the argument used when v was an eigenvector 


we can show that D(v, u') = 0. Consequently 
D(v, t) — »,H(v, t) = 0 


where t is an arbitrary admissible vector. This may be written as 


(i d't! (; iy) ( 1s .) ( bi fs) 
i (i x <;)-B (i x ds? ds + B et xX ds? “= * B oe | x ds? : 


al 


— Xr | mA(i X v)-(i X t) ds = 0. 
Integration twice by parts yields 


. 2 d°v ) 
[ ‘5 Ee x Fy — mAX.(i X v)?-(i X t) ds 
us J 


, ( 1 = da” : ( ~ ? qd? 
+ (i x i). (i x rs) | - E x t): “ (B ix ) | = 0. 


Since the vector t is chosen arbitrarily in 0 < s < l and also at s = I, it follows im- 


mediately that 


d’ . E i es 
We E ix ad — mAy,i X v = 0, 
as s as 


and ats = l, 


piix@io 4x a = 0. 
ds ds \ ds 
These are identical with the differential equation (2.10) and the boundary conditions 
(2.12). Furthermore \; = A3-1 2 *** 2 Xi , for the class of functions used in minimizing 
D(v) for 2 is more restrictive than any previous class. Also it can be readily shown 
that \2 is actually the nth eigenvalue in order of increasing magnitude. For if uw is an 
\2_, and t is the corresponding eigenvector then 


eigenvalue different from eg Kes 
“ee D(t) 
e H(t) 
and H(t, u’) = Ofori = 1, --- n — 1. Thus t is an admissible vector in the minimizing 
problem considered. Therefore 
D(t) > Div) 
H(t) ~ H(v) 


or 
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Consequently A; is the nth eigenvalue. 
Finally D(v) and H(v) may be expressed in scalar form as 


al 
D(v) = | (EL,Q® + EI,P?) ds + [e(E1,Q)* + €(EIsP)"|, , (3.7) 
al 
H(v) = | mA(vi + v3) ds, (3.8) 
where 
,_ ar, dv, dr 
silica ds* er ds ds V2 7 
(3.9) 
a d’v, 9 dv, dr — 
o-"a* "ata" e™ 


The differential equation (2.8) and the energy principle derived in this section can be 
extended to include torsional vibrations about the i direction coupled with transverse 
vibrations. If the axis of twist and the axis of centroids do not coincide we obtain coupled 
transverse-torsional oscillations. If they do coincide, the differential equations of motion 
uncouple and we have torsional oscillations superimposed on the transverse vibrations. 

4. Rotating twisted beam. The vector differential equation derived in Section 2 for 
transverse vibrations will now be extended to a twisted beam which is rotating. The 
bending moment M,(s) due to the internal forces, which might be called the elastic 
moment, must be equal to the moment due to the external forces. This moment is the 
sum of three moments; M,(s), the moment due to the inertial forces relative to a rotating 
frame; M,(s), the moment due to rotation; and M-(s), the moment due to the Coriolis 
forces. 

In order to compute these moments the right handed coordinate system Xy , X, , X2, 
with unit vectors I, J, K is introduced (Fig. 2). The I axis coincides with the i axis. The 


XJ XJ 
ie J 


Q 
LOIN Q KL XK 


LO) 9) 


Ko) Ke A 
x, A 


SO) 





Fia. 2 


beam is assumed to rotate with constant angular velocity 2 about the J axis. If 6 is 
measured as indicated we have 








252 R. C. DiPRIMA AND G. H. HANDELMAN [Vol. XII, No. 3 


G(s) = [ rds+ 60. 


“0 


When the beam vibrates the centroid of the cross section Q moves to a new point 
Q’. The position vector R(s) of Q’ may be written as 


R = sI + U,(s)J + U,.(s)K, 
= si oe u,(s)j + uo(s)k, 
where the transformations from one coordinate system to the other are 


U, = u, cos 6 — usin 8 u, = U, cos 6+ U2 sin 6 


U, = u, sin 6 + uz cos 6 Us 


II 


—U, sin 6+ U, cos @ 


(4.1) 
J = cos 6j — sin 6k j = cos 6J + sin 0K | 
K = sin 6j + cos 6k k = —sin @J + cosOK ) 
If a is the total acceleration of the centroid Q, it is easy to show [3] that 
aU, aU. 
a, = Sy + oe 
ot ot 
ar = —O{sI + U.(s)K} 
ac = 20 ou, 
at 
Multiplying a by — mA gives the effective force intensity at the cross section, and 


hence the moment at the cross section s due to the action of an infinitesimal element 
of the beam at &, — > s, is 


—mA([R(t) — R(s)] dé X a. 
If second order terms are neglected, that is, products of the U’s and their derivatives, 
the I component of M;(s) and the Coriolis moment, M,(s), must be neglected [4]. Con- 
sequently we set 


M,(s) = M,(s) + M,,(s) (4.2) 
where 
M,(s) = mA(é)(E — 8) .. £ aoe dij — / mA(é)(E — 8) 8 (6) dé K 
J ; at” 
ai 
M-(s) = 2? | mA(é){t(U.(€) — U.(s)] — (€ — 8)U.(E)} dt J (4.3) 


~@ | mA@eHU@ — UW) a K.| 


“s8 
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From Section 2 we have 


M,-(s) = Bi y4 ar 9 


Ss 


Q 
cS 


tol 


Q 


Noticing that dJ/ds = dK/ds = 0, and differentiating Eq. (4.2) twice with respect to 
s we obtain 
P du . a1 aU, ; 
2. (m4 ¢—) = of | a Li | — mA(s)l sohy 
s \ IS os 

(4.4) 
a'Us(8) 


at K, 


+ 2 : E 10) |x + mA{s) A. 2(8) J — mA(s) 
Os Os ot 
where 


al 


L(s) = | mA (g)é dé. 


Either the capital U’s and J, K can be transformed into the small u’s and j and k or 
vice versa. Doing the former we obtain as our equation of motion 


3 J 7°u)| . : Ju 
.{Bo-i x <5? — o fix! 
Os | Os 


As Os ) 
: (4.5) 
— mA(s)%@-(i X u) + mA(s)i X u = O, 
where the dyadic @(s) is defined as 
cos’ 6 —sin 6 cos 6 
O(s) = 
—sin 0 cos 6 sin’ @ 
If the beam is cantilevered about both axes at s = 0 the constraining boundary 
conditions are simply 
: —— 
ixu = 0, x as = 0; (4.6) 


and, since the moment and its derivative due to the rotation vanish at s = l, the re- 
quirement that there be no shear and no moment transmitted across the free end requires 


er 8 less SR 
B *s x as’ —_ 0, as {B 1 x ay _— 0. (4.7) 


Assuming a solution of the form 
u(s, t) = v(s)e 
(v,(s)j + v.(s)k]e™ 
= [Vi(s)J + V.(s)K]e™ 


and substituting 


ws) = i XY, 
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we obtain the following differential equation 


oy ; : , 
# {g.4 wl - 9 é Lr) at — mA(s)2@-w — mA(s)’w = 0. (4.8) 


ds’ | ds") 1s 


The corresponding scalar differential equations for the displacements v, and v, are 


@ (EI,P) — 2r ¢ (ET,Q) - qt ET Q) — 7°(EI.P 
os “" ds ah® — i a all 


: lv, a , 
_ of 119 PC) — mata) (4 _ r,)s — v, sin” #@ — v, sin 6 cos |} — mA(s)rX'v, = 0, 
as 


@ EI Q) + 27 @ ny p) + 2 (BLP) — EI Q) 
ast | “" dg ©"? } ds LL of) T (Ld, 


. lv, P ' , 
~ of 1900 -~ marco ( + 7) — v, sin @ cos 6 — v2 cos” a} — mA(s)rX\v, = 0, 
ds 


where P and Q are defined by Eq. (3.9). 
The boundary conditions are 


l 
aia a = (4.9) 
at s = 0, and ats = / 
2 j re 
B ee = 0, ¢ (B is 4 = 0. (4.10) 
ds ds ds 


The variational principle stated in Section 3 may be generalized to give a minimum 
principle which yields the eigenvalues and corresponding eigenfunctions of Eqs. (4.8), 
(4.9), and (4.10). The functional D’’(w) = D’’(i X v) which plays the role of the potential 


energy is 


ve wf dv aia ’ 
D”"Ud X y) = | (i x Y). IM z(S) — M’7(s)} ds, (4.11) 
“0 as 
where 
M’i(s) = B-i x 2S 
E ds”? 
and 
M'k(s) = 2 | mA(é&){é|V.(&) — V.(s)] — (& — s)V.(E)} dé J 


—- ¥ | mA(é)E[V(E) — V,(s)] dé K 


al 


= —2L(s)i X v(s) + & | mA(§)(i XK v) dé 


—- 2 mA(E(E — 8) V2(s) dé J. 
S/\S : 


vs 
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Substituting these expressions in Eq. (4.11) and integrating by parts in the proper 


manner we obtain 


i fiw _ dw 


pw- | a+ SS 
J, ds 


: ! 
ih mA()Yw-O-w/ ds, (4.12) 


ds° 


or in scalar form 
Di X v) = | {ET,Q” + EI,P? + x10) (2 + ro) 4b (ae -- w.) | 
0 s s 
— mA(s)(v, sin 6 + v, cos a} ds. 


The functional H’’(w) which plays the role of the kinetic energy is 


H’"(w) = [ mA(s)w-w ds (4.13) 
or in sealar form 
H(i X v) = [ mA(s)(vi + v3) ds. 
We introduce the mixed forms 
Dw, t) = [ 4 4 + YL(s) aw at - mA@w-o-t} ds,  (4.12a) 
and 
Hw, t) = [ mA(sw:-t ds (4.13a) 


Since @ is a symmetric dyadic both D’(w, t) and H’’(W, t) are symmetric functionals. 
It also follows immediately that 


D’''(w + t) = D’'(w) + 2D’, t) + D’(t), (4.12b) 
H''(w + t) = H’'(w) + 2H’’(w, t) + A’(t). (4.13b) 


The vector w(s) will be called an admissible vector if it is continuous and has con- 
tinuous derivatives of at least the fourth order, and also satisfies the constraining boundary 
conditions Eq. (4.9). It must be pointed out that due to the presence of the minus sign 
in the last term of D’’(w) it is not obvious that D’’(w) is positive for all admissible w, 
and it appears to be difficult to prove that such is the case. A careful investigation of 
D"'(w) reveals that the troublesome term arises from V,(s). It is felt, but not justified, 
that this term is always at least balanced by the other terms in the functional. 

The general minimum principle is the following: Of all the admissible vectors w(s), 
that one which minimizes D(w) under the side conditions H(w) = 1, and H(w), t') = 0 
fori = 1, 2, --- , nm — 1 where t’ is the 7th eigenvector of the differential equation is an 
eigenvector of Eq. (4.8) and satisfies Eqs. (4.10) as natural boundary conditions. The 
minimum value of D(w) under these side conditions on w is then the corresponding 
eigenvalue \2 and \2 = AZ_, = +--+ = Aj. The proof of this principle is analogous to the 
proof of the variational principle given in Section 3, and will not be repeated. 
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An orthogonality relationship between the eigenvectors of Eq. (4.8) can be developed 
in the usual way. If w(s) and t(s) are two eigenvectors, and \” and y’ their respective 
eigenvalues where \° ~ yu’, then w and t satisfy the orthogonality condition 


H’'(w, t) = | mA(s)w-t ds = 0, (4.14) 


which implies, as before, that 
D’'(w, t) = 0. (4.15) 
The reality of the eigenvectors and eigenvalues of Eq. (4.8) can now be proved in a 


manner analogous to the proof given in Section 3. 


If we assume rt = 0, and set v. = 0, H” and D” reduce to 
t Sd Lv lv, dv 
a at 7 f . dv, a . ee 
D'(@,) = | (+ (21, ) + 2L(s) _: mAQ'v; sin” @> ds, 
: q ds * ds ds ds J 


H'(v,) = | mAv; ds. 


Applying our variational principle we obtain the differential equation 


d’ f ‘ d’v, | La dv, : af . 
3 \ EI, 3" — & Lis) 7 — mA», sin” 6. — mAXv, = 0, 
ds? \"*? ds? ds ds | 
and 
dv d {,,, dv 
== (EI, ) = 0 
ds ds ds 
at s = |. These are the differential equation and boundary conditions developed by Lo 


and Renbarger for a rotating blade vibrating flexurally in a plane making an angle of 
(x/2 — 6) with the plane of rotation [5]. 
Finally we wish to develop a similarity principle which may be used for experimental 


work. Introducing the dimensionless variables 
z= s/I, y = w/1, 6 = cl. (4.16) 


Eqs. (4.8), (4.9) and (4.10) become 


/ 2 2. i 
| ) . ‘B™ f : i x) S} — mAy-y — mA ay = QO, 


07 1*7 dx’ | ‘da*{ ~ dx ( 
with 
== () dy 0 
y 7 , dx — 

at xz = 0, and 

1 bea a 

p< = 0 (B o2) = 0 | 

dx dx ax 

at x = 1. Also 


. 


- cos y —sin y cos 
K(2) = | mA(n)n dn; y= 
—sin y cos p sin yp 
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with 


Vix) = | ddn+ 0. 
“0 
Consider another beam which we designate by C’ which has the same @ and 6 as 
the first beam, but which is rotating with an angular velocity 2’. If the dimensions of 
C’ are related to those of C by a similiarity factor a, then the area A’, and moments of 
inertia J/ and J; of C’ are related to those of C by A’ = a°A, I{ = a‘l, , I; = a‘I, , and 
if l’ = el, E’ = bE and m’ = cm, then we have for the beam C’ 


‘ba? d’ ei d’y d ‘ ay\ — 
ew F x a _ - —=—>-?7 — Aw: a A = = . 
(a 7) dx {B ty dx {kK (x) dx mA€-y mA Q”” atti 


The homogeneous boundary conditions remain the same; consequently if 


bat _1 
ce OQ’? ~ 9”? 
then 
a »* 
~ 


5. Numerical example. As an example, the natural frequency of a twisted, uniform, 
non-rotating, cantilevered beam will be computed using the energy principles developed 
in Section 3. For simplicity we assume the beam has a constant natural twist, that is 
7 equals constant. The energy principle may be written for v an admissible vector, 


’ s D(v)/H(y). (5.1) 


Writing D(v) and H(v) in scalar form, introducing dimensionless variables, and using 
the fact that the beam is uniform we obtain 





2 mM Al’ 2 
dion EI, ‘. 
ta (a? Ww, dw. 2 \ ’ (eu dw, 2 \ (5 .2) 
ge! ow ’ —-,;- 2 a> eum , 
- | dz? 26 ya bw,? dx +y [ qe? + 26 i 5'w,? dx 


al 
[wt + wh ae 


where y = J,/J, . It is known that 6’ = 12.36 for a beam with no twist [6]. 

Of particular interest is the calculation of the natural frequency of a propeller blade, 
in which case the ratio of width to thickness of the beam is very large, that is, J,/J, is 
large. Use of the variational principle in this case requires some delicacy in the choice 
of the approximating functions. Since y is large, h, = (i X 8°u/ds2), , which is called the 
first bend, will be small. Consequently we shall choose trial functions which make the 
contribution to the potential energy from the h, term small. Now 
du 


— “a 2 _ 
M=B 1X 3 
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and expanding this we obtain 
M, = EI,h,, 
s Ou» Ou Or : 
EL 5 - = oF + — th, tied 
Os 5 s 


From this we see that we wish to choose trial functions that will make the coefficient 
of y in Eq. (6.2) as small as possible. Let us choose w, and w, as fourth degree polynomials, 


w (2) = xX + ast + AyxX , 
. 4 ; 2 
w(xr) = fu + dx + ba’. 


For the particular computations in question, the constants a, and a3 were determined by 


satisfying the conditions 


aw . d ‘d?w, , 
—- Ow, = 0, (‘ =~ — sw,) = 0 


dx dx \ dx” 

at x = 1 for 6 = 7/16. These equations result from requiring that there be no moment 
and shear in the j direction upon assuming w, = 0. In this case we obtain aj = —3.974054, 
a, = 5.980107. The constants b and d are determined from the conditions that the linear 
and constant term in (d°w,/dx° + 2édw,/dx — 6°w,)° should vanish. This gives 

9 

p= § d= —= éa,. 
Silat 


The constant f is left arbitrary to be used as a minimizing parameter. If p is defined as 
the ratio of the natural frequency of the twisted beam to the natural frequency of the 


untwisted beam, we obtain the following results* 


6 I,/T: 6? p 
x/16 48 12.65 1.01 
64 12.72 1.01 
144 13.07 1.03 
r/8 48 13.39 1.04 
64 13.69 1.05 
144 15.01 1.10 
1/2 48 14.21 1.07 
64 14.64 1.09 
144 16.85 1.17 


These results appear to agree closely with experimental results [7,8]. To obtain better 
approximations for larger values of 6 and /,/J, , w, and w, should be chosen to be poly- 
nomials of higher degree, thus allowing us to require that the coefficient of y in Eq. 
(6.2) vanish to a higher degree in z. 


*The authors wish to express their gratitude to Mr. Howard Montgomery of Carnegie Institute of 


Technology for his assistance in carrying out the computations. 
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Some of the computational aspects of this problem are being investigated by Henry 
E. Fettis of the Flight Research Laboratory, Wright Field. He has noted the delicacy 
of the calculation problem and has proposed transformations of the energy principle 
which appear to improve the situation. He has also investigated several different com- 
putational schemes. These results will be published at a later date. 
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ON THE AXIALLY-SYMMETRIC STEADY WAVE PROPAGATION IN 
ELASTIC CIRCULAR RODS* 


BY 
JULIAN ADEM 
Instituto de Geofisica, Universidad Nacional de México and 
Instituto Nacional de la Investigacién Cientifica, México, D.F. (Mexico) 


Summary. In the first part of this report we find an exact solution for the problem 
of steady wave propagation in an isotropic, elastic circular bar of infinite length, free 
of stress on its lateral surface and loaded by a harmonic body force, parallel to its axis, 
whose amplitude is a Dirac delta function. The solution of this problem gives at the 
same time the solution for a semi-infinite bar with a special prescribed load on its bound- 
ary plane. 

In the second part, we find a solution for the semi-infinite bar in which the conditions 
at the boundary plane are prescribed in terms of functions which give implicitly the 
stresses and displacements. 

Finally, using a frequency of interest in current ultrasonic experimental work, we 
develop a numerical example and compare the result with the case of a low frequency. 

1. Infinite bar with the body force 6(z)e"‘*‘. We consider an infinite circular bar 
of a perfectly elastic isotropic material, free of stresses on the lateral surface and loaded 
by the body force 6(z)e~**‘, where 6(z) is the Dirac delta function defined by 


d(z) = 0 for z #0, 


[ (2) ds = 1, 


and w is a positive number. 

The problem is that of determining the displacement at all points of the rod for the 
steady case (i.e. the time dependence for stresses and displacements is e~***). 

The general solution. For an elastic isotropic medium the equation of motion is 


(A + u)V(V-u) + wV*u + pX = pu”, (1) 


where u is the displacement vector; X is the body force vector, \ and yu elastic constants 
and p the density. 





» (DB 











Fic. 1 


Using cylindrical coordinates (r, 6, z) as shown in Fig. 1, we have axial symmetry 
(i.e. the solution is independent of @). 





*Received Aug. 11, 1953; revised manuscript received Dec. 21, 1953. The results presented in this 
paper were obtained in the course of research sponsored by the Watertown Arsenal under Contract 
DA-19-020-ORD-2598 when the author was at Brown University, Providence, R. I. 
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Letu = Vo+ V7 X Fand X = ki(z)e"**', where k is the unit vector in the direction 
of the z axis and V-F = 0. 
Assuming that 


or, 2, t) = ofr, ze '*’, 
F(r, z, 2) = Fr, 2e *”' 


and considering that 
V[S(z) + C] = ké(z), 


where S(z) is the unit step function defined as 
S(z) = 0 for 2 = @, 
S(z) = 1 for z2> 0, 


and C is an arbitrary constant, we have that (1) is satisfied if 


> 


Ve sis hg = —([S(z) + C] a (2) 
V(VX+k(V XF =O, (3) 
where 
42 = —Pw_ 2 pw 
r+ Qu’ . 


The formula dS(z)/dz = 6(z) has been used in the past in a formal way, without rigorous 
justification. However, recently its validity has been proved on the basis of the theory 


of Distributions. The reader is referred to L. Schwartz, ‘‘Théorie des Distributions’’. 
Using the formula V?(V X F) = —V X V X (VY &X F) and assuming that the 


only component of F that does not vanish is that in the @ direction (which will hence- 
forth be called F), Eq. (3) becomes: 
9 [1 00rF wF . 2 
‘ k O(rk | 4 df 4 iF = 0, (4) 
or Lr or Oz 


while (2) written in cylindrical coordinates becomes (taking C = 0): 


\ 


1 0 dy 0” 9 S(z)h° - 
—{r——)4+—$4+ ho = —-— >. (5) 
r or or Oz @ 
S(z) can be expressed as 
’ 1 12 e'*? 
Si) = = | —— de. 
271 ¥ g 


where the path of integration is the real axis, indented to pass below the origin. 
Let 


l ite ‘ 
g(r, 2) = on ; e*" e*(r, ) dé, 
(6) 
: aa 
Fir, z) = on / e'**F*(r, £) a, | 
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where the path of integration must be such that the solution represents outgoing waves 


at infinity. 
From (4) and (5) we obtain: 


Ly eae 22") 2 32) mh i a 

2r | ‘ E or (; or + Ch ee" + itw” pe = 0, 
Lp” [21208 4 a — ere fo ae = 

2r | | 2 r or + (h oF | & = 0. 


These equations are satisfied if 


h- 





1d 3e") - en. 2 
r or ( or ray = ea 
| 7) 
7 9 (J Ki 
91 ArF*) | pepe — 0, 
orr or 
where ao =? —$,B8 =k’ -—& 
A solution of these equations is: 
h? | 
¢* = = 2 oo A(é) J o(ar), 
Lew a ( 
| (8) 
F* = BE)J,(8r), 


where A() and B(~) must be determined from the boundary conditions. 
Without loss of generality we assume that the radius of the rod is 1; therefore the 


boundary conditions are r,, = 0,0, = 0, atr = 1 
But 
1 O(ru,) , du Ou, 
a : = the ~~ oe 
7 NE or + as +> or’ 





Ou, Ou, 
Tre = Al de + a 


and the displacements in the r and z directions are 


(9) 


Therefore our boundary conditions become 


; Wo - 2 o°F 
{ata 340% 4085-25) =o, 


Fare) 0 2 : 
9 - = § = 
{2 Or dz ~ Oz aiel r=1 
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and in terms of ¢* and F*, 











0°o* a P _ OF* 
lo + Qu) =a +4 . AE g* — gu a ] = 0, 
or or Oe 34 
(10) 
2 , de" + (2¢? = ky ‘| 
“1g 07 <—§S . .= 0 
Substituting (8) in (10) 
A(&)2itad (a) + Blé)(k? — 2)J,(8) = 0, 
; : NEN? 
A(§) {Quad ,(a) — [(AX + Qua” + AE ]J.(a)} + BE)[—2iuBEJo(B) + 2tud,(B)] = —¥4. 
wa 
The solution of this system of equations is 
Nhe? J (a) 
Bg) = -—— 
() 2uw a 
| (11) 
; ) ER me 2 ETL 
A(t) = _ ah (3k ~—— ES (8) | 


2uw aA 


where A = £a8J;(a)Jo(8) — a3k?Ji(a)Ji(8) + (3k? — £)?J.(a)J,(8). Substituting 
(11) into (8), and the resultant values of (8) into (6), we obtain 





Se Da MO inh? f® &(—#? + 4k’) Ji(8) J (arye* .. 
fil alias ~ Drie” | . tar =- 4 3 in a A is, 
(12) 
2 p@o 2 tke 
P(r, 2) = — DEL f° ELI 
4uw 7 J ad 


In order to compute the integrals that appear in this solution it is convenient to deter- 
mine the roots of A/S = 0. 

The frequency equation. The equation 
JB) _ 


2 
= faJd,(a)Jo(8) — = : kK? J (a) J,(8) + (2 ke — ?) J o(a) 8 0, (13) 


a 

B B 

where a = (h? — #’)'”, B = (k*? — 2)”, known as Pochhammer’s frequency equation, 

has been the subject of several papers*, because it gives the modes of longitudinal wave 

propagation in an infinite circular rod, having stress free surface (velocity of the wave = 
w/t; wave length = 27/). 

In this equation, given the elastic constants, ¢ is a multiple-valued function of the 
angular frequency w. For the purpose of computing our integrals we are interested in 
determining all the roots of Eq. (13), corresponding to the given frequency w. 

Let us consider the frequency equation (13) for large | ¢|. In this case a and 6 can be 
written as 

a = it — h*/28), 
B > i(€ — k*/28) 


*D. Bancroft, Phys. Rev., 59, 588 (1941); A. H. Holden, Bell System Technical J., 30, No. 4 (1951); 
R. M. Davies, Trans. Roy. Soc., (London), A, No. 821, 240 (1948). See also the references mentioned in 





these papers. 
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and the frequency equation for sufficiently large | ¢ | becomes 
2i¢ — cosh 2§ = 0. (14) 
Let & = & + 7t, , where &, and &, are real; then Eq. (14) can be written as the pair 
2, + cos 2% cosh 2¢, = 0, 
—2¢, + sin 2& sinh 2¢, = 0. 
Considering that | € | is large, these equations have solution only if | £, | is large and, 
since 
£,/sinh ¢, > 0 when |t, | @, 


it can easily be shown that there is an infinite number of roots of Eq. (14) which are 
approximately at the points 


& = +3 log 2kx — ikr/2, 
£ = +} log 2px + tpr/2, 


where k is a large positive even integer and p is a large positive odd integer. 
In Eq. (13) the variable is £’; therefore we must consider the solution for £ as well as 
for —£. For the minus sign case a and 8 are given by 


: h? : k? 
a> —ie(1 ~ ft), = ~ie(1 _ #) 
and the frequency equation becomes* 
2i¢ + cosh 2§ = 0, 


whose roots are those of (14) with opposite sign. 
Therefore the frequency equation (13) for large | ¢ | can be written as 


(27 — cosh 2£)(27¢ + cosh 2¢) = 0, 
or 


42” + cosh’? 2¢ = 0 


and there are an infinite number of roots of it which are approximately at the points 
1 _ 
£ = +5 log 2kr + ix, 


where k is a large positive integer and where all combinations of sign are allowed. 
Since there are no roots of (13) on the axes for large | ¢ | it follows that there are a 


finite number of roots on them. 
Finally, we are going to show that there are no roots in the neighborhood of the 


real axis (except at the axis itself). Let us denote (13) by the brief notation 
Fé) = 0; (13’) 
this equation can be written as 


oii F(,) + ifeF'G) + +--+ = 0 (13"”) 


*As was also shown by C. W. Curtis in a private letter to the author. 
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and, if £, is small enough to neglect £3 , since F’(£,) * 0 when F(é,) = 0, and F(é,) and 
F’(é,) are both real, it follows that (13’’) can be satisfied only if & = 0. 
Evaluation of the integrals. From now on we shall consider the solution for z > 0. 
The integral 


can easily be evaluated by using Cauchy’s residue theorem.* There are three poles on 


the real axis: one at £ = 0, and the other two at = +h; in the rest of the complex 
plane the integrand is regular. To satisfy the condition at infinity (outgoing waves) 
we choose the path of integration as shown in Fig. 2, i.e. from & = —R to& = R indented 
at = +h and é = 0 to pass above the real axis for § = —hA and below it for the other 
poles; and then from ¢ = R to = —R along a semicircle I with center in the origin. 


Integrating on the chosen path and letting R — ©, since (by Jordan’s lemma) the 
integral on the semicircle [ goes to zero we obtain: 


= 7 [2 —e™). (15) 





Now let us consider the integral 


] ee Fr e&(-& + tk’) J (8) J (are acl dt: 
2 a A Ss? 


“ ) 


dividing the numerator and denominator of the integrand by 8 we obtain: 


I, = | f(® J olar)e - dé, 


where 


-_— L 2 J,(8) 
a-t +35*}", 


hee <p 


The numerator and the denominator of the integrand are both regular functions and 
the only singularities of the integrand in the finite plane are simple poles corresponding 
to the roots of the equation a’A/8 = 0. 


*See E. T. Whittaker and G. N. Watson, ‘‘A course of modern analysis’’, Chapter VI; or any of the 


standard books on complex variable. 


1954 WAVE PROPAGATION IN ELASTIC RODS 267 


Let & , f35. °° , be the positive real roots and &,,, , £42, -** the complex roots of 
A/sp = ( 

We shall again use Cauchy’s residue theorem, choosing the path of integration as 
shown in Fig. 3. 


4S, 


4 | 

’ A Tr ! 

' | 

} | 

} | 

‘ a ae | 

een en eer ae ny! $1 

Fia. 3 


The path goes from A to B on the real axis, indenting it at the poles to pass above 
the real axis for negative poles and below the axis for positive ones (the reason for this 
choice of the path is to satisfy the condition of “outgoing waves” at infinity); then 
we close the contour of integration by a family of paths [, such that 


lim | f(®) J o(arje'* dt = 0. 
Elmo /Tp 
In our case we have used the family of rectangles shown in Fig. 3: the side CE passes 
between two consecutive pairs of roots of the frequency equation (for large | & |) and 
CD = DE = EA = 3x(N + 3), where N is a large positive integer. 
Integrating on the chosen path and letting N — © we obtain: 


| f()Jolarje"®* dt = 2ni > R, — lim f(@) J olar)e’® dé, 

J —@ q=1 N-o /Tp 

where R, are the residues at the poles on the positive real axis and upper half-plane. 
Sut, since 

lim | f(®J fare dé = 


Noo / Tp 


evaluating the residues we obtain: 


[ f() J olarve : dé = 2m bis (Eq) J o(agre'*** 2. we ihe ae 
where 
1p.2 _ ¢2 
g(t) = SOE Ae) 
dz (4/8) 


and a,g= [a], 6 
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In a similar way the last integral of (12) can be evaluated, and we obtain: 


[7 ELLIE ge ons F uaanuoere™, - 
where 
i(t) = E Fifa) _ 
aB dé (4/8) 
and 
Bo = [Blenee - 


Substituting (15), (16) and (17) in (12) we see that the terms in e’* cancel and obtain: 


1 rh? 


g(r, 2) = —3+5—3 Dd g&.)Jo(arde*", 
w 2UW on 
| (18) 
. . Nh7i . - / ik oz 
a we 2 UE) Ji(Barje*". | 
Substituting (18) in (9) and considering the time dependence, we obtain 
u(r, 2, t) = 4 > {&,GE,)Jo(ar) — B.T(E.)Jo(B.7) je" *”, (19) 
q=1 
u(r, z,t) = >> {—a,G.)Ji(ar) — &.T(E)J(8.7) eS", (20) 
q=1 
where 
2 £/ _ ¢* 2 .” ; 
G() = =, : + Mate | 
PANTO 2, a - 
a B dé (A B) 
(21) 
vn . We #Ji(a) 
T®) =5-3- pte , | 
Zw ; ; 
af dé (A/B) 
&, , &, °°: , & are the roots on the positive real axis of A/B = 0} E41 , &s2, °°* are 
the roots in the upper half plane; 
a= (h? = ¢”)! ¥ 8 a (ke? = . 
1/2 1/2 
pe(ey™, 2- (2) 
A + Qu my 
a = la}gee, » B, = [B]e-e, - 


In the solution (19), (20) the terms due to the complex roots have the decay factor 
e~** (“a” real and positive) and for z sufficiently large, they become very small and can 
be neglected. As we consider the solution for smaller z, some of those terms probably 
must be considered (those that correspond to the poles that are closest to the real axis). 
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We will obtain better information about the characteristics of the solution (19), 
(20) in the last section where we discuss a numerical example. 

To obtain the solution for z < 0 we integrate as in Fig. 3, following the same path 
AB, but instead of the path BCDEA we take its reflection in the real axis. It can easily 
be checked that u,(r, z) = u,(r, —2), u,(r, 2) = —u,(r, —2z), as is to be expected by 
considerations of symmetry. 

2. Semi-infinite bar. Consider the same bar as in the previous problem. The equation 
(1) without body force is 


(\ + pw) T7(V-u) + uV7u = pu”. (22) 


We look for a steady axially-symmetric solution of Eq. (22), which satisfies the 
stress-free condition at the lateral surface of the bar and at the plane z = 0 is of the 
type u, = u,(r)e**', u, = u,(r)e*** [or o, = o,(r)e***, te = Tra(rle **'). 

The solution (19), (20) of the above problem suggests that we assume: 

u=Vet+VvV XF; V-F=0, (23) 
where 
io, 


F(r, Zz; t) = F(r, aati 


I 


g(r, 2, t) 


where £ is to be determined from boundary conditions. Substituting (23) in (22) we 
obtain 


1a se) om 
12 (,2¢ ayes, 
d 1 A(rF) > 

Orr or ver +S 


where F is the @ component of the vector F, the other components being zero; and a 
and 8 have the same meaning as before. Therefore 


g(r, —) = AJ (ar), 
F(r, §) = B&)J,(6r). 
From the boundary conditions at r = 1 
1 ’ a 2 


ttaJ (a) ’ 


where ¢ must be a root of the frequency equation (13). 
From (9) we obtain 


+ 





u, = . Bs Bee.) [—§.M(E,) Jo(a@.r) + Bod o(Ber) Je’ “E"~*”, 
(24) 


u, = § D BE)—aME (ar) — BG """”, | 
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where 


M€,) = 
ts £,a,J (a,) 


and é, , a, , 8, have the same meaning as in solution (19), (20). 

Similarly we can express o, and r,, in terms of the arbitrary functions B(é,). 

In order to satisfy the conditions at infinity (outgoing waves of finite amplitude), 
in the solution (24) we must include only the terms corresponding to the roots of the 
frequency equation that are on the positive real axis or above the real axis. 

At the plane z = 0 we can prescribe displacements or stresses by prescribing the 
functions B(é,). To do this we need to know the roots, £, , of the frequency equation (13). 
If we know only some of the roots, then for the other roots we can take B(é,) = 0. 
We can in this way construct solutions in which at z = 0, u, and u, are given by (24); 
(or o, and r,, by the corresponding formulas). 

The problem of an infinite bar with the body force 6(z)e"**’, where 4(z) is the Dirac 
delta function, corresponds to the case in which B(é;) = —7T(é;), where T(é;) is given 
by (21). 

3. Numerical example. Let us consider a bar of radius equal to 1 cm. 

We shall take Poisson’s ratio equal to 1/3: vy = 1/3. For this case* V,; = 2V, , where 
V, = (A + 2u/p)’” and V, = (u/p)'” are the velocities of the longitudinal (compres- 
sional) and transversal (shear) waves in an infinite medium. 

Therefore k = 2h. 

Let 


E = hy, 25) 
then the frequency equation (13) becomes: 


y(l — y)'? fh — y/)'* J ofh(4 — y')'7] 


2(1 — y°)”* ThA 2) "1d, [h(4 pyres 
— wae a Y. 2 — 1 : ’ ~ 
h4 — y*)” il Yy il Y) 
+ Boek oD 2\1/2 2\1/2 ; 
+ (4 mis Jolh(l — y/) “| J, [h(4 — x)" = O. (26) 
(4-—y 


For the moment we are interested in the case of frequencies, w/27, in the ultrasonic 
range which vary approximately from 5 Xx 10° cycles/sec. to 25 X 10° cycles/sec. To 
be specific we are going to consider the case w = 27 X 10’; furthermore we shall take 


9 \1/2 , 
v. os (2 + 2s) = § X 10° cm/sec., 
» 4 


then 





4 P Gold. 3rown Universitv Techn teport W A-6487 /4. 
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The positive real roots of the frequency equation (26) are in the interval 0 < y < 2.16. 
There are approximately 130 positive real roots. Those in the interval .85 < y < 2.16 
are given in Table 1. 


TABLE 1 
Roots of the frequency equation (26) for h = 40x in interval: 0.85 > y > 2.16. 








8514 1.0006 1.6414 1.9263 
. 8609 1.0336 1.6592 1.9332 
8715 1.0733 1.6753 1.9397 
8835 1.1090 1.6914 1.9459 
.8949 1. 1437 1.707 1.9517 
. 9046 1.1775 1.7222 1.9572 
9126 1.2119 1.7368 1.9624 
9208 1.2418 1.7509 1.9672 
9291 1.2725 1.7646 1.9717 
9387 1.3018 1.7778 1.9758 
9466 1.3257 1.7905 1.9796 
9534 1.3571 1.8030 1.9843 
9589 1.3846 1.8147 1.9862 
9645 1.4099 1.8262 1.9893 
.9708 1.4344 1.8373 1.9915 
9759 1.4583 1.8479 1.9937 
9810 1.4817 1.8581 1.9956 
9855 1.5038 1.8680 1.9971 
9894 1.5251 1.8774 1.9983 
9926 1.5464 1.8865 1.9992 
.9954 1.5664 1.8952 1.9997 
.9974 1.5863 1.9035 2.1500 
.9989 1.6056 1.9115 

.9998 1.6236 1.9191 





The roots on the imaginary axis closest to the real axis are at 
y = .1591, 


y = 192i. 


The roots &, are 407 times the corresponding roots of (26). Substituting them in 
the solution (19-20), we see that of the terms corresponding to real roots the ones with 
largest amplitudes are those that correspond to the roots near y = 1. If in the solution 
we consider only the terms corresponding to the roots in .86 < y < 1.4, we get a result 
with an accuracy roughly speaking, of 1%. 

The terms corresponding to the imaginary roots have the decay factor 


—40r\uglz 
e q 
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where the smallest value of | y, | is 0.159, and therefore such terms can be neglected in 


the solution, except in the neighborhood of z = 0. If we take z > 1/2, for example, we 


can ignore such terms. 
The solution, considering only the roots in the interval .86 < y < 144, is 


u,(r,z, t) =i a > [TJolagr) + S.J o(B,r) Jerr?” 


> I,(a,r) Y f - jn 
* » |e, I,(a,) + SAAB) ¥ ; 


— 4 _ 1 f 7 ,1(40 ry qgs—wt) 
u,(r,z, i) = Or? p> M J ,(B.r)e 


where 7/, , @ » Be » T'¢; Se» Yo s %p » By» Rp, Sp, M, are given in Table 2. 


TABLE 2 


( 


7) 














Qq Be le Sq M, 
. 8609 63.9127 226.8361 + .02 +.05 
8715 61.7387 226 . 2329 — .09 + .28 
8835 58.7479 225.4412 +.11 +.15 
8949 56.0587 224.7626 —.10 +-.17 
9046 53.4575 224.1217 + .05 +.17 
9126 51.8238 223.7448 + .02 — .22 
9209 48.9587 223.0913 —.11 — .16 
9291 16. 5082 222.5761 +.14 —.18 
9387 43.5677 221.9729 —.15 — .16 
9466 10.7403 221.4325 +.13 —.18 
9534 38.0762 220 . 9676 — .08 —.13 
9584 36.0279 220 . 6032 — .03 + .28 
9645 33.4141 220.2010 +.15 +.31 
9708 31.039 219.8617 —.27 + .35 
9759 27 .3571 219.3716 +.21 + .22 
9810 24.3914 219.0198 — .23 +.23 
9855 21.4760 218.7056 + .25 + .24 
9894 18.2590 218.4166 — .27 + .30 
9926 15.2305 218.1778 +.30 +.38 .20 
9954 12.3151 217.9893 — .32 + .46 26 
9974 9, 2363 217.8385 —.45 + .80 44 
9989 6.1575 217.7380 — .45 + .94 54 
9998 5.0266 217.6752 +.77 +2.06 1.18 
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I a, Bp R, Sp M, 
1.0006 1.3480 217.6123 + .42 —7.31 —4.07 
1.0336 34.2937 214.9357 + .08 +1.55 92 
1.0733 18.8205 212.1083 + .06 —1.20 —.73 
1.1090 60.2308 209. 1552 + .03 +.79 54 
1. 1437 69.7184 206 . 1769 0 — 61 — 36 
1.1775 78.1253 203 . 1484 0 + .47 36 
1.2119 86.4945 199.7304 0 — .21 —.18 
1.2418 92.5641 196.9909 0 +.31 
1.2725 98.7718 193.9490 0 — .29 
1.3018 104.7786 191.1978 0 + .22 
1.3257 110.2827 187.6415 —.11 
1.3571 115.2716 184.6255 +.19 
1.3846 120.5997 181.1949 —.11 
1.4099 125.0859 178.1162 +.11 
1.4344 120.3334 175.0500 —.11 
1.4583 133.3295 172.034 +.11 
1.4817 137.6146 168.616 —.11 
1.5038 141.1706 165.6629 +.11 
1.5251 145.0163 162.3076 — .07 


Some remarks. The longitudinal displacement for z not near the origin (say, z > 1/2) 
consists of approximately 130 harmonic waves (one for each positive real root of the 
frequency equation) which travel with velocities varying from V,/.118 to V,/2.15 
and with wave length varying from (1/20)(1/0.118) to (1/20)(1/2.15), respectively. Since 
in the solution (27), Jo(a,r), Jo(8,r) and Io(a,r)/Io(a,) have absolute value less than 1, 
the coefficients T, , S, , R, , S, , M, give the order of magnitude of the amplitude of the 
different waves. The waves of largest amplitude are those which travel with a velocity 
approximately equal to V,; = [(A + 2u)/p]'” and with wave length approximately 1/20; 
the closer the velocity is to V, the larger the amplitude of the wave is; and it is interesting 
to point out that the amplitude of the wave that travels with the velocity V,/1.0006 is 
remarkably larger than the others. Therefore if we are looking for an approximate 
solution, it is sufficient to consider only the roots near the point — = h(i.e. y = 1). 

An example with low frequency. Now let us consider the solution for small values 
of w. 

When w — 0, h — 0, k — O and the frequency equation (26) has only one root at 
y 3/2 and therefore the solution will have only one term. 


The wave velocity 


and the wave length tends to infinity. 
Let us consider h = 3 (w = 3V,). In this case the frequency equation has 3 positive 
real roots which are at y = .692, y = 1.080, y = 2.105 and the solution (not near the 


origin) {considering the imaginary part of (19), (20) only] is: 
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u, = 4V7° [.397J,(2.1651r) + .574Jo(5.6286r)] cos (2.0762 — 3V,2) 
+ 4V;* [—.004/,(1.8792r) — 0.604./,(3.8043r)] cos (3.5402 — 3V;t) 
+ 4V_,~ [—.008/,(5.5569r) + 0.0117,(1.9695r)] cos (6.3152 — 3V,0), 
u, = 4V;" [—0.406/,(2.1657r) + 0.021./,(5.6286r)] sin (2.0762 — 3V,2) 


LV >* [—0.0027,(1.8792r) — 0.562./,(3.8043r)] sin (3.5402 — 3V,2) 


LV>? [—0.0071,(5.5669r) + 0.086/,(1.9695r)] sin (6.3152 — 3V;,t). 
In this solution the amplitude of the two first terms decreases as r increases, while that 
of the last one increases and gives practically the only contribution near the surface of 
the rod. The velocity with which the wave corresponding to the last term travels is 
almost that of Rayleigh surface wave. 
The wave lengths are 


2r 2r 2r 
L,=——, lL, =z, L, =>. 
2.076 3.940 6.340 
The wave velocities are 
V , V, : V 
eee, Keg, Fimo. 
0.692 1.1 2.105 
Figure 4 shows the displacement vector at several points of the section z = 5 cm 


for h = 40x and for h = 3, for ¢ an integer. 

General conclusions. The solution (24) [or (19), (20)] consists of an infinite number 
of terms of which all but a finite number (those corresponding to the positive real roots 
of the frequency equation) decay with z. 

In some of the following statements we shall consider » = 1/3, but the conclusions 
apply to the general case. 

When w (frequency) increases, the solution is changed as follows: 

1) The number of waves that do not decay with z increases. For example if in the 
formula w = AV, we take h equal to 1, 2, 3, 3.75, --- We obtain 1, 2, 3, 4, --- waves, 
respectively. For h = 407 (ultrasonic frequency) we obtain 130 waves. 

2) The rate at which the terms corresponding to the complex roots decay with z 


increases. 

3) Since the wave lengths are given by 27/hy,; , where y; are the positive roots of 
Eq. (26) which are in the interval 0 < y; < 2.5, as w increases the wave length decreases. 
In our example for h = 407, the wave length of the terms of larger amplitude is approxi- 
mately 1/20, while for h = 3 it is 27/2.076. 

4) The roots of the frequency equation (23) become closer together particularly at 
&=h(y = 1) andé k (y = 2), and the roots near these two points correspond to 
waves that travel with velocities that approach V,; and V, respectively as w increases. 
The velocity of the wave corresponding to & > k (y > 2) approaches the velocity of 
ayleigh surface wave. In our example for h = 407, the terms of larger amplitudes are 
those with V = V, 

5) At a section z = constant, the displacement vector will vary in amplitude, direc- 
tion and sense as r varies. Such variations become larger and faster as w increases. 


(As shown in Fig. 4). 
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The amplitude of the wave corresponding to the root §; > k decreases when r de- 
creases, and it is practically a surface wave. The amplitudes of the other waves decrease 


as V increases. 
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ON THE NON-STEADY MOTIONS OF A RIGID BODY IN AN IDEAL FLUID* 


BY 
G. W. MORGAN 


Brown University 


1. Introduction and summary. [If a rigid body moves in vacuum under the action of 
a single external force which has a constant direction but possibly varying magnitude 
and whose line of action passes through the center of mass of the body, the body will 
maintain a pure translatory motion parallel to the direction of the force. If the body 
is surrounded by a fluid, (in the following assumed to be inviscid and incompressible), 
such a motion will, in general, no longer be possible. The following question arises: 
Is there a point such that a body acted on by a force which has constant direction and 
whose line of action passes through this point, will move in a straight line parallel to the 
force without rotation? The investigation presented below gives the following answers 
to this question. 


a) There is no such point in the most general case. 

b) If the shape of the body satisfies certain conditions, then there exist three axes, 
each associated with the body and with one of three directions relative to the 
body, (the so-called directions of permanent translation)**, as well as with the 
density of the fluid, such that, if a force of constant or varying magnitude acts 
on the body with its line of action along one of the axes and if there are no other 
external forces, the body will move with pure translation parallel to the force 
and axis. The conditions which the body shape must satisfy are tantamount to 
certain symmetry requirements. 

c) In the case of plane motion (without circulation), i.e. if the body is an infinite 
cylinder with its generators perpendicular to its direction of motion, there exists 
a unique point (for a fixed fluid density p) such that a body, whose motion is 
due entirely to an external force which acts parallel to one of the two directions 
of permanent translation with its line of action passing through the point, will 
move without rotation parallel to the force. This point might be referred to as 
the ‘apparent center of mass” of the body corresponding to the given fluid. 

d) For a body of revolution there also exists an apparent center of mass. 


Inasmuch as these considerations concern a problem in classical hydrodynamics, 
one might expect to find them in the standard literature on the subject, but, to his 
surprise, the author was unable to do so. 

Our plan of procedure is to set up the expressions for the kinetic energy of the body 
and the fluid and then to investigate the equations of motion of the body in terms of 
the kinetic energies. 

2. Kinetic energy of the body. Consider a system of rectangular Cartesian axes fixed 
to the body with origin at O. Using tensor notation, the motion of the body can be 
described in terms of the translational velocity U, of the origin and rotation with angular 

*Received Aug. 19, 1953. This paper represents a slightly condensed version of a report prepared 
under Contract N7onr-35810 between the Office of Naval Research and Brown University. 

**Lamb, Hydrodynamics Ch. VI. 
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velocity w; about the origin. If x; denotes the position vector of a point of the body 
with respect to O, then the kinetic energy of the body, 7’, , is given by 
oT, = | (U5 4A €5550;4,) (Us + €:1,@:2,)0 dr, (1) 

/B 
where o is the density of the body, e¢;;, is the alternating tensor, and the volume inte- 
gration extends over the entire body B. Since 

€cjx€cinUeLn = (6; :%,X, — X;X1), 
this can be written as 


Zi, = UU, | o dr + Ujye,eij, | ot, At + w,U <€;;; ot, dt + ww,l,, , (2) 
“B ~“B 


“B 
where /;; is the inertia tensor of the body and is given by 


Il, = | o(6;;X,.X, — x;x;) dr. (3) 


/B 


We note that this expression involves ten independent constants: 


the mass of the body | odr=M, 

“B (4) 
the three moments | ot, dr = Mz} , 

“B 


where x* denotes the position of the center of mass, and the six independent components 
of the symmetric inertia tensor /;; . 
At this point it is convenient to introduce a new notation. Let 


V; = OP , V2 = U2 ; V; = U;, = 
(2) 


V, = Wy V; —= We; V6 @3 . 


By V, we shall mean any of the six components V, to V, with @ ranging from 1 to 6, 
We use Greek subscripts to differentiate the present notation from tensor notation. 
Also, supposing the summations implied in the expression (2) to be carried out and the 
U, and w; to be replaced by the appropriate V, , let M., be the coefficient of the term 
in V,V, . The 36 quantities M,, then have the following values: 

M,, = Mz = M;; = M, 


Mi, Mo; = M;; - Ms, = Moz; —_ M32 = 0, 
Mis = May, = Mos = Moz = Mae = Moz = O, 


M,, = Ms, = —My = —My. = M27, 

M,. = M,, = —My = —M, = —Mot, 
(6) 

M., = Mya = —M;, = —M;, = Mz, 

M 4 _ Is; ’ M,, = Fes 9 Mee = | 2 

Pe Mas = Toy; hence M,, = M;,, 


Mus = Ts, Me = Ig, hence My, = Ma, , 
M.. } Bax Isc ; hence M,, = M,; . 
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We can now write expression (2) in the compact form 
6 
2T , _ : B VaV 5M as . (7) 
a,B=1 


3. Kinetic energy of the fluid. Since the motion of the fluid is assumed to be entirely 
due to that of the solid, it is irrotational. Its velocity potential is of the form 


ee TVs, (8) 
a=l1 


where the ¢, are functions of the coordinates x, , x, , x; depending solely on the shape 
and position of the solid surface relative to the coordinate axes fixed in the body. They 
represent the velocity potentials due to the motion of the body with unit translational 
or rotational velocity with respect to the appropriate axes.* The kinetic energy of the 


fluid, 7, , is given by 
- D i  - 
27, =p | &— ds, (9) 


where p is the fluid density, 0/dn denotes differentiation with respect to the outwardly 
directed normal to the surface S, and the integration extends over the surface S of the 
body. 

Substitution of (8) into (9) gives 


oT; = > (V.V sp | ¢ Ovs as). (10) 


ares on 


Each integral 
F O08 16 
ro, — aS 
I - on 
is a constant determined exclusively by the shape and position of the body surface S 
relative to the coordinate axes. Hence, setting 


Ap 
o | va dS = Fas, (11) 
Js on 


we can write (10) in the form 


6 
oT = >» VaV ek ap . (12) 
a,Bp=1 
Since V’o, = 0 for all a, 

F Ogg Y [ Oa Y 

Da dS = re ae dS 
I, re an Js ¥8 “an 
so that F., = Fs. , and hence Eq. (12) involves 21 independent constants. 


4. Equations of motion of the body. The total kinetic energy T of fluid and body 


combined is 


T=T,+Tr =} Dd, VaVpAes; (13) 


ap 


*Lamb, loc. cit. 
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«ec 


where A,g = M.s + Fag . The equations of motion of the body in terms of T' are: 
d eT ae. oe 
dtav, — Vo as Vs oy, =X, 
d eT oF . OT _ 
dtaV. Veoy, + Vo OV, = X,, 
d oT re i . oT 
dtav = ‘oy. * ‘av. - Xs, 
{ | +) 
d oT - OF. mee gy ee oc : oT _ 
dtav, — ‘av. V2 5V. ~ Vay, 7 Vs ay. = 4s, 
d oT . ar ee oo a OE * 
dt av. - 4 ev « T Ot , *% ‘OV, + Ve OV, = Xs, 
d oT _— . oe 0 . or “ 
dt av, a“ av, y , 0 V. * : OV, ¥ : OV, = X ms 


where X, , X, , X; are the components of the external forces and X, , X; , X, are the 
components of the external moments with respect to 0 that are applied to the body. 
(By external forces and moments we mean forces and moments other than those due to 


the fluid pressures. ) 

5. Rectilinear motion. We now choose a special orientation for our coordinate axes, 
namely that of the three mutually perpendicular axes of “permanent translation’’ of 
the body. These have the property that if the body is set in motion parallel to one of 
them, without rotation and in the absence of external forces, it will continue to move 
in this manner; i.e., for motion with constant velocity in one of these directions the fluid 
exerts no force or moment on the body. We shall now proceed to examine the external 
forces and moments necessary for accelerated motion of the body in one of these direc- 
tions. 

It can be shown that if the directions of permanent translation are chosen as co- 
ordinate directions, the constants 


F Fi, = FF, = F. = Fy = Fs. = 0. (15) 


In the following we shall always consider this choice to have been made. 


6. Motion with velocity V,(¢). Let us consider first the case when V, = V,(¢) and all 


other V,, = 0 for all time. Then the equations of motion (14) together with the expression 
for the kinetic energy (13) give: 
V,A,, = X VA, = X, VA, X V,Ay = X, 
(16) 
ViA,; — V;A Y; ViAie + ViAi \ 
where the dot denotes differentiation with respect to time. Hence, from (6) and (15) 
X Se. CAA Le X, = 0, 
17) 
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Expressions (17) represent the system of forces X, , X. , X; whose line of action 
passes through 0, and the system of moments X, , X; , X, about 0, which must be applied 
to the body so that it may move in the prescribed manner. Our object is to ascertain if 
there exists a point 0°” such that, if the line of action of the forces X, , X, , X3 passes 
through it instead of through 0, the external moments X$'’, X;'’, X,'’ with respect to 
0°" will vanish. If such a point exists, then, under the sole action of the external force 
X, applied through 0", (we may now conveniently call this force X;"’), the body will 
always move without rotation and in the direction of the x, axis. 

To investigate this question we refer all quantities to a new coordinate system whose 


1 


origin is at a point 0°’ and whose axes are parallel to those of the original system. 


ae 5 A . ‘ r(1) 1 1 
Let the velocities of the body referred to this new point be V{’; and let d;"’, d2", 


1 


d;,’ be the components of the displacement of 0°’ from 0. Then 


$ = r. Vs = ie Ve = emg 
V, e VEO st Vode = VX”, 

(18) 
V, = Vs? + VePds? — Verdt”, 


V,= V;’ + Vi"d; — Vad’. 


If we substitute expressions (18) into the formula (13) for the kinetic energy and 
collect coefficients of VS'’V3", which we eall A} , we obtain 


Aw - Aj, + Ajod3"’ = A,;3d;"’, 
Ais = Ai; + A,3d;"” — A,,d,"’, (19) 
Axe = Aig + A,d;" -_ Awi,’. 


The remaining coefficients are not of interest at the moment. 

If we now rewrite the equations of motion with respect to the V{'’ and 0“? we obtain 
equations identical with (17) except that all quantities will now have superscripts (1). 
Hence the external moments about 0°’ can vanish only if 


Ais _ Ais 


= Ai, = 0. 

Using (19), (15) and (6) this means we must have: 
F,, = Fi,’ = 0, 
Fis + Mis — (Fir + Mids” = 0, (20) 
Fig + Mig + (Fix. + M,,)d2” = 0. 


Since F,, does not vanish in general, equations (20) cannot generally be satisfied 
and hence there exists no point with respect to which the moment of fluid pressures 
vanishes for accelerated motion in the direction of the x, axis. It is clear that analogous 
conclusions will be reached for motion along the two other axes. 

We may note that this result could have been deduced directly from Eqs. (17) 
referring to the original system. The components of the external force which can con- 
tribute to the moment X{” are X, and X, ; but these vanish. Hence we must have 
X, = X;"’, i.e., the moment is the same for any reference point and is therefore due to 
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a pure couple. Thus a change in the point of application of the external force can have 
no effect on X,. 

Returning to Eqs. (20), we see that if /,, happens to vanish for a particular body, 
then all Eqs. (20) can be satisfied by choosing 0‘ appropriately. Using (6) we then have: 


ny , 1 ~* ee Coe 1 2 
dq = Pis + Mz’ BP iis Fis Mx? (21) 
F,, = M F, + M 
Hence in this case we can find a point 0“, or rather an axis 2,"’, (since d,' is arbi- 


trary), such that under the sole action of a force X,') along this axis, the body will 
always move parallel to this axis and without rotation. We might call such an axis, 
an “axis of rectilinear translation’. 

The condition that F,, = 0 imposes a restriction on the body shape and so might be 
regarded as a symmetry requirement. 

If the original origin 0 is taken at the center of mass of the body then the terms 
Mzx* and Mzx* vanish in (21). 

7. Motion with velocity V,(t). Repeating the above considerations for the case when 
all V, except V. vanish for all time by transforming to a new system with origin at 
0”, we find that the external moments about 0 will vanish only if 


Fox + Mog + (Foo + M22)d3” = 0, 
F, = Q, (22 


Fog + Mog — (Foo + M22)d;”’ = 0. 


As expected, these equations have no solution in the general case since F,, is not 
generally zero. 

If F,,; happens to be zero, Eqs. (22) are satisfied if 
Fy, — Mz3 q? = Fr6 + Mri | (23) 
F.,+M ’ F,. + M 

Hence there exists an axis 7."’ analogous to the axis z,'’ found previously. We note, 
however, that the two axes will not, in general, intersect since d;'’ ¥ d;” except under 


a a —_ 


special circumstances. 
Since an analogous result will be obtained for motion with velocity V,; we can sum 
up our results as follows: 
(i) In general, rectilinear translational motion under the sole action of one external 
force is not possible. 
(ii) If the body shape satisfies a certain symmetry condition, e.g., F,, = 0 when 
referred to coordinate axes parallel to the axes of permanent translation, then 
there exists an axis of rectilinear translation for motion with velocity V, 
Similarly, if F.,; = 0 or F3, = 0, there exist axes of rectilinear translation for 
motion with velocity V, or V; , respectively. 
(iii) In general no two of these axes, if more than one exists, will intersect. 


8. Plane motion without circulation. If the body is an infinite cylinder with its 
generators parallel to the x; axis which moves with velocity components V, , V, and 
V, only, then the fluid motion will take place in planes parallel to the x, , x. plane and 
the two axes of permanent translation will lie in this plane. As before we choose our 
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coordinate system in the direction of these axes. For such plane motion all coefficients 
A.» in the expression (13) for the kinetic energy having one or both subscripts equal to 
3, 4 or 5, vanish. 

Considering motion with velocity V, only, we obtain equations (17) with 
A,, = A,; = Oand hence X, = X; = 0. The transformation to a new coordinate system 
with origin 0"” yields one equation for d;'’ which will make A{}’ = 0. It is the same as 
the second Eq. (21). This determines one axis of rectilinear translation. 

Similar arguments hold for motion with velocity V,(t) and they lead to the second 
Eq. (23) for d;” which will make A;;’ = 0. This determines the second axis of rectilinear 
motion. 

Since both axes lie in the plane of motion, they intersect at a point P, say, with 
coordinates x;"’ and x,"’. If the origin of our first coordinate system is taken at the center 
of mass of the body, then x;” and x," are given by the equations for d\”’ and d}"’, respec- 
tively, with z} = z$ = 0: 

(a) Foe (a) Fie 
2, = F.+M’ 2 = F.+M’ (24) 
We might call this point the “apparent center of mass of the body’’, keeping in mind, 
however, that its position also depends on the density of the fluid which appears in the 
constants Fag . 

It is interesting to consider pure rotation about this point with velocity V,(¢). The 

equations of motion are found to be: 


A\? Ve — ViAS? = Xi, As? Ve + VaAi? = X:", 
(25) 
X,=X,=X,=0, ASV, = Xe, 

where the superscript (a) means that we are referring all quantities to a coordinate 
system whose axes are the axes of rectilinear motion and whose origin is therefore z;°, 
x,"’. Since A;s’ = A,{’ = 0 the only external force or moment on the body is a moment 
about the origin. Thus the’point 2x;*’, x,” also acts as an apparent center of mass in the 
case of pure rotation about this point. 

It is important to note that the motions and corresponding forces and moments 
considered above cannot be simply superimposed. A translation of the apparent center 
of mass, for example, with velocity components V, and V, but without rotation, cannot 
occur without the action of an external moment about the apparent center of mass. 
This is due to the non-linearity of the problem. The analogy between the apparent 
center of mass for motion of a body in a fluid and the center of mass is restricted to three 
particular motions, the two translations of the body parallel to the axes of permanent 
(or rectilinear) translation, and rotation about the apparent center of mass. 

9. Body of revolution. The axis of revolution is a direction of permanent translation. 
Let it be the x, axis. For motion in the 2, direction we consider Eqs. (20). For a body 
of revolution F,, must always be zero because changing the sign of V, in the expression 
for the kinetic energy must leave the latter unchanged. Similarly F,; = F,. = 0; also 
M,; = Mxt = 0 and M,, = —Mzx* = 0 since the center of mass lies on the axis of 
revolution; hence ds’? = d;'? = 0. This gives the obvious result that the axis of rotation 
is an axis of rectilinear translation. 

From symmetry considerations it is clear that the choice of directions for the axes 
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x, and x; is immaterial. For motion with velocity V, , say, we must consider Eqs. (22). 
Applying symmetry considerations to the expression for 7', as was done above, we find 
that F., = F,, = 0. Hence 


>) 


Ags = Ma + (Fo + M2)d3” = 0 


(26) 
Ase = Fo, + Mo = (Fx. + M,.)d;"’ = 0. 
But 
M. = —Mzx* = 0. 
Thus d;”’ = 0, and there remains 
* Fi, + Mzx* = 
i" «= , (27) 
F.2 + M 
which defines an axis of rectilinear translation. 
For motion with velocity V,; we obtain in the same manner d;” = 0 and 
7 F,, — Mx* 
di ae Z as > (1 “ (28) 
F33 + M 
From the symmetry of the body it is seen that F,. = F;, and F., = —F;,; . Hence 


as was to be expected from our assertion that the choice of the directions 


x, and x; was immaterial. 

Thus, as in the case of plane motion, there exists a center of apparent mass, but 
here there are infinitely many axes of rectilinear motion: the axis of revolution and any 
axis perpendicular to the latter and passing through the apparent center of mass. 

Let us consider a motion of pure rotation about an axis of rectilinear motion. Clearly, 
for motion with V,(¢), i.e., rotation about the axis of revolution, the only force or moment 


e _ d? 


is 
X = MY V,. 


For motion with velocity V, , say, about the appropriate axis of rectilinear motion we 


have: 
( a r2 (a) r(a) 
Ais Ve — VeAcs = Xi", 
Aze Ve + VeAis = X2°, 
Azs Vo = X;°, 
(29) 
(a) 77° 72 (a) ry (a) 
Ais < gee VeAse _ X4 ’ 
Ass Ve + VeAce = X5°, 
( "ad r (a) 
Age 6 - Xe . 
But Ajs? = Ags’ = 0 and A3%’, Ass’, Ass’ vanish from symmetry considerations. Hence 


all forces and moments with the exception of X;” vanish. Thus the concept of the ap- 
parent center of mass for a body of revolution is also valid for pure rotation about an 
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axis of rectilinear motion in general. An exception to this is the motion due to a com- 


bination of velocities V,(¢) and V,(é). 
As in the plane case, however, simple superposition of the special motions and forces 


discussed is not permissible. 

It should be noted that, even if all assumptions made in this analysis are sufficiently 
well justified, some of the special motions discussed cannot be expected to be observed 
in nature because they are unstable when subjected to small disturbances. 
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ON THE AVERAGE SECOND MOMENT OF THE ENERGY SPECTRAL 
INTENSITY* 


BY 
TSE-SUN CHOW 


Scientific Laboratory, Ford Motor Company 


Introduction. In a paper by J. Leray [1] it was proved that the kinetic energy of a 
two-dimensional incompressible viscous flow in a domain D bounded by a regular curve 
S dissipates at least as fast as an exponential law. Following this, J. Kampé de Feériet [2] 
proved the same proposition by using a different technique and related the dissipation 
constant with the domain 9. In a later paper [3] he studied the Fourier transform of the 
two-dimensional vorticity equation and proved that the upper bounds of the Fourier 
transform of the vorticity and the energy spectral function also decrease in accordance 
with a similar law. 

The present paper is concerned with an incompressible viscous fluid such that the 
domain of flow extends to infinity. One of the objects of this paper is to study the time 
variation of the average second moment of the energy spectral intensity which is closely 
related to the dissipation of energy. 

We begin by assuming the velocity and pressure field to belong to a class of functions 
such that they can be represented by a trigonometric series. For this purpose we introduce 


K(r-x) = cos (r-x) + sin (r-x), (1)** 
so that 
v(x, 2) = >> A(t, )K(r-x), (2) 
and 
p(x, t) = >> &(r, dK(r-x), (3) 


where v(x, ¢), p(x, #) are the velocity and pressure respectively, being functions of the 
position vector x in the physical space and the time ¢, A(r, ¢) is the velocity spectral func- 
tion and is a vector, ®(r, ¢) is the pressure spectral function, and ris a wave number vector. 
The purpose of using K(r-x) as defined by (1) instead of the customary exponential 
function exp (r-x) is to avoid complex spectral functions in (2) and (38). 

The fundamental equations governing the motion of an incompressible viscous 
fluid are the Navier-Stokes and the continuity equations: 


= +v:-Vv=- : Vp + Vv, (4) 
ot p 


V-v =0, (5) 


where V is the gradient operator and p, v are constants, the density and kinematic vis- 
cosity of the fluid. Using (2) and (3), we can transform (4) and (5) into 


< A(r, ) = S[A(r, )] — AG, d, (6) 


*Received Sept. 30, 1953. 
**This notation, K(r-x) = cos (r-x) + sin (r-x), was first introduced, to the author’s knowledge, by 


L. S. G. Kovasznay. ris a vector of magnitude r. 
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where 
‘ » owe r ' a 
F[A(r, t)] = 4] QO, ) — 3 QC, O-r (7) 
| r 
and 
Qir, ) = >> [Air +s, ) + A(t —s, #) + A(—r+s, t) — A(—r—s, #)]-sA(—s, ft). (8) 


Af, t)-r = 0. (9) 


A study of the transformed Navier-Stokes and continuity equations (6) and (9) has 
been made by the author [4]. Here we summarize a few of the conclusions which will be 
made use of in this paper: 

(a) Let the spectral function A(r, t) be defined by the spectral vectors A, , A, ,A;,... 
.in the wave number space at time ¢, then S[A(r, ¢)] represents the 


located atr, ,f.,T 
dt at time ¢ due to the non-linear interaction of the spectral 


contribution to dA(r, ¢) 
vectors A, , A, A; ,... ; S[A(r, ¢)] is in general not zero at time ¢ forr = +r; +1, , 
i~j,i,j = 1,2,3,...andiszeroforr + +r; +1; . It is to be noticed that S[A,(r, ¢) + 
A,(r, t)] + S{A,(r, t)] + S{A.(r, £)] in general, while S[cA(r, )] = c’S[A(r, t)], c being a 
constant. 
(b) For an arbitrary spectral function A(r, ¢) which satisfies the orthogonality relation 
(9), A(r, t)-r = 0, the following is true: 
> A(r, t)-S[AC, t)] = 0. (10) 
r 
(c) >> A(x, t)-A(r, t) = E(t), where E(é) is twice the average kinetic energy of the 


flow per unit volume. In what follows, we shall call E(¢) simply the energy of the flow and 


A(r, ¢)-A(r, ¢) the energy spectral intensity. Therefore, by combining (6) and (10), we 


obtain 
ad. re 
7 E(t) = —2v >> r’ACr, t)-A(r, d), 
1.e. 
E(t) = E(O) exp ee | (r°(1) ar | (11) 
where 


> rA(r, t)-A(r, t) “ae 
(12) 
>> A(r, t)-A(r, d) 


(r(t)) = 


and is the average second moment of the energy spectral intensity. 

From (11) it is seen that the average second moment of the energy spectral intensity, 
{r’(t)), is intimately related to the energy of the flow, E(t). The paper is mainly devoted 
to the study of the time variation of this average second moment. 


1. From (6) we have 


~ A(r, t)-A(t, t) = 2A(r, )-S[A(e, 0] — 2°AC, 0)-AG, t). (13) 
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It is seen that the energy spectral intensity A(r, ¢)-A(r, ) varies with the time ¢ due to the 
interaction of spectral vectors and due to viscosity. It follows that the average second 
moment of the energy spectral intensity may also vary with ¢ due to the interaction and 
viscosity. We shall, first of all, examine the effect of viscosity alone, ignoring the effect of 
the interaction. The following conclusion can be easily deduced: 

I. The average second moment of the energy spectral intensity is a decreasing function of 
time if only the effect of viscosity is taken into consideration and will approach r, as t + 


where To = min (fr, , Te, 13 5++>)s 
When the interaction term is ignored, Eq. (6) becomes simply 


a ° 
; / At, ) = —w°AG, 0), (14) 
Cc 
so that 
A(r, t) = A(r, 0) exp (—r°d). (15) 
A set of spectral vectors initially located at r, , r. , 13 ,... Will remain at r, , 1% ,13,... 


at any later time, each behaving individually and damped exponentially. Thus 


: > °° A(r, 0)-A(r, 0) exp (—2vr*t) 


r | t) = - : 2 (16) 
p A(r, 0)-A(r, 0) exp (—2pr 0) 
and it can be verified that 
d | 
aa* 
2. A(r, , 0)-A(r; , O)A(r; , 0)-A(r; , O)[r? — 77]? exp [—200% + rE] 
Dy — — LS (17) 


21 a ere 
{>> A(r, 0)-A(r, 0) exp (—2" 0)} 
Statement I is thus seen to be true. 

2. When viscosity is disregarded the energy of the flow is conserved, as is obvious 

from (11). Then from (12) we have 
( 2 f | d 

(Tr (t)) = 

dt > A(r, t)-A(r, é) dt 


Thus it is sufficient to examine d/dt bm r’A(r, t)-A(r, 0). 
It follows from (6) that 


> Ar, t)-A(r, 2). (18) 


l . . 
= > r’A(r, )-A(r, ) = 2 >> Ar, d)-S[ACr, DO], (19) 

( 
from which we have: 

II. The time rate of change of the second moment of the energy spectral intensity due to the 
effect of the interaction of spectral vectors alone is zero if at that instant the spectral vectors 
are so located atr, ,t,%,,...that tr; +r; +r, #0, forall possible i, j,k = 1, 2,3,..., 
LF | x k, k #14. 

It is to be noticed, however, that although in this case the first derivative of (r°(¢)) 
is zero at this instant, the second derivative is not zero in general, due to the fact that 
the interaction of spectral vectors starts to generate new vectors at +r; + Tf; . 
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We proceed to examine the more general case for which the velocity spectral vectors 
are so located that one or more of +r; + r; + r, may be zero. Consider, for the simple 
case, that there is only one such set, say, r, +r. +1; = 0, then from (7) and (8), 


J 


F[A(r, t)] = Jy [A(r +s, t) + Afr —s, ) + A(—r-+s, d — A(-—r-—s, #)]-sA(—s, 0) 


S > [A +s, 2) + Afr —s, 2) + A(-—r +s, #) — A(-r —s, #)]-sA(—s, ry, 


r 
so that 


r’A(r, t)- S[A(r, 2)] 


Z zz {fA +s, 2) + Afr —s, 2) + A(-r-+s, t) — A(-r —s, #)]-sA(—s, 2) 
-A(r, 0)} 


r >. ‘Ar —r;,0 + Ar +r,,t) + A(—r-—r,, — A(-—r+r,;, d] 


-r,A(r; , 2)-A(r, td}. 


Writing A(r; , 2) = A; , we see that 
TAs TA, ‘A, + rsA.°T,A;-A, 


D) 


> A(t, t)-s[AC, )] = —3 + r3A,-r.A.:A; + riA;°TA,-A.;. (20) 


+ riA.-T3;A;°A, + r2A, ‘T;A,-A, 


We evaluate the right hand side of (20) by choosing the coordinate axes so that 


r, , rz all lie in the plane of R,-R, , as shown in Fig. 1. The angie between r; , r; is 








I, , 
Ry 
as 
A; 
Y; 
car © 
5 ; | ya Rs 
| - | JA 
mei i f 
| ~ 
R A; 
Fic. 1. 


denoted by 6;;, measured from r; to r,; in the clockwise direction when facing the R;-axis. 
bof ; .* ‘ 
Since the velocity spectral vector A, is perpendicular to the wave number vector r, 
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because of continuity, its direction is completely fixed by the angle a; which is measured 
from a parallel vector to the R,-axis to the spectral vector A; in the counter-clockwise 
direction when facing in the direction of the corresponding wave number vector. Using 
this choice of axes it can be verified that: 


(sin (82, — 63) COS a; COS a2 SiN az 





J ° ° 
) + sin (0:2 — 603) COS a; SIN a COS a3/. 


) rA(r, t)-S[A(r, )] = 3A, AcAsrirers | 
a | (21) 


ri.? 





+ sin (63; — 6:2) SiN a; COS a» COS a; } 
Equation (21) is also true for +r; +r; + r, = 0. Obviously, we can generalize the result 
so as to include the case in which there are more sets such that +r; + r; + r, = 0, 
i,j,k = 1,2,3,---,t1# 9,7 #k, k # 1. Hence we have, 

III. When only the effect of the interaction of spectral vectors is considered, the time rate 
of change of the average second moment of the energy spectral intensity is given by the following 


formula: 


er 
r(t)) = . 
sin (6;, — 6;) COS a; COS a; SIN a, | 
Zz A,A;Arrarin ) + sin (6;; — 6;,) cos a; sin a; cos Oe (22) 
+ sin (6; — 0;;) sin a; COS a; COS ay, 


where the summation is extended to all sets of vectors r; , 1; ,¥, such that +r; +r; +r, = 0, 
i,j,k = 1,2,3, °° ,t #797 Ak, k #2. 

Nothing can be said about the sign of d(r*(t))/dt for the general three-dimensional 
flow. However, if the motion of the fluid is two-dimensional, then d(r*(t))/dt is zero and 
r(t)) remains constant with time. Hence: 

[V. The average second moment of the energy spectral intensity of an unbounded two- 
dimensional incompressible non-viscous flow remains constant with time. 

3. We now examine the time rate of change of (r’(t)) due to the combined effect of the 
interaction of spectral vectors and the viscosity. As a matter of fact it follows from previ- 


ous results that 


— ae . ° 
isin (6;, — @,) COS a; COS a; SIN a, 
d J . . 
r(t : pm A,;A;Ajrit jn + sin (0;; — 0;,) COS a; SIN a; COS a, | 
dt a Yee... 
| + sin (6,; — 6;;) sin a; COS a; COS a, 


ia AT LA, -A;A;-A,(ri — 7°)’. (23) 


Therefore we have 


V. The average second moment of the energy spectral intensity of an unbounded two-di- 
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mensional incompressible viscous flow is a decreasing function of time; therefore, there exists 
a lowe r hound of energy: 


E(t) = E(O) exp { —2v(r°(0))é}, (24) 


where (r°(0)) denotes the initial average second moment. 

For the three-dimensional flow it has not been possible to establish a similar energy 
bound valid for all time ¢t. However, the following shows that sueh a bound exists at 
least for the initial stage of dissipation. Starting with (6) we have 


1 ‘ 
= > r°A(r, t)-ACr, 2) + 2 > r‘ACr, O-AC, 0) 
= 2 >> 7°ACr, 1)-S[A(r, 0] 


< 2{>> r‘A(r, t)-A(r, ) >> S[ACr, ]-S[A(r, d]}!. 


Hence 
| 4 , j F ( i ( , | 4 
4 & Pag, )-AG, ) + ATE MAG, 0-AG, O} | oN Bie o1 , 
dt \ 4y° J 
- dL» SIAC, #)]-5[AC, 0] 
= Dy , 
and therefore 
d vos F[A(r, t)]-F[A(r, 2) 
FT pm r A(r, t)-A(r, t) S 2» Sl 0] bd — J . (25) 
“av 
Here we need a lemma 
LemMA. Let A(r;,f) = A;,7 = 1,2,3,..., then 
> slA(r, )]-s[Ar, 0] $2 4 A? rr A?. (26) 


Proof. It has been shown by the author [4] that given two spectral vectors such that 


A(r, ,?) = A, , A(t. , t) = A, then 


: (ri — 73)" sin’ a, sin” a, 
A, Az | sin O12 (g — —poge 





| clas 4 7 1 
is A te I = 3 Bs 2 ¢ 
? Ti + 72 + 2rir2 COS O12 
‘ y | 
+ (r,; COS a; SIN G2 — 1, COS a SIN a)” (27) 
forr =r, +r,andr = —r,—-f, 
. , (ri — ci sin’ ay, sin A 
S[A(r, t)] | = $A,Az2 | sin 62 | {2 7 > 
iri + 72 — 27,72 Cos 8,2 
; 3 
+ (r; COS a, SIN ay + 72 COS a2 SiN a)” (28) 
forr =r, — fr, andr = —r, + 1,, where the wave number vectors r, , r, lie in the plane 


of R,-R, , and the various angles are measured using the same convention as before, Fig. 2. 
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R3 





,-% af | 





Using these results, we have 
> s[A(r, t)]-F[A(r, 0] 
< 14?A4% sin’ 6,.[2(r, + 72)? sin’ a, sin? a. + 2r3 sin? a, cos’ a, + 2rj cos’ a, sin’ ag] 
= A}AX(r, + 12)° 
S 2A; Ar; + 72). 
Now let A(r; , t) = A;,7 = 1, 2,3,..., then 
> s{A(r, )]-5[AG, ] S$ 2 D> AAri + 7) 


37.412. 49, 


IIA 


and the lemma is proved. 
From (25) and (26) we have 


> r’A(r, t)-A(r, t) S ye r’A(r, 0)-A(r, 0) exp f | b A(r, 7)-A(r. 7) ar}, (29) 
Jo 
1.€. 
1 t 
r°(t) < r°(0) exp {2 | E(r) ar}, (30) 
0 
where r’(¢) is the second moment of the energy spectral intensity. We have, then 
@ (1) => —2y°(0) exp {t [ E(2) ar} (31) 
dt os o : ; V Jo ' 


from which we have 


VI. For an unbounded incompressible viscous flow the energy is bounded below given by 


E(t) = BO)! — 2/” oo | exp {EO ‘ _ i} (32) 
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and the average second moment of energy spectral intensity is bounded above given by 


» (r°(0) | {E(0) ,\ im {Z0) \ 
r(t)) < ((0))\1 — 27? > xp {—— t? — xp 4—— 33 
) Ss t v E(0) exp | m ty | exp ¥ ty (33) 


for 


Il 


y E(O) | 
< a ; 3 
;' E(0) in [ + 2v*(r° (0) | (34) 
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—NOTES— 


ERRORS IN ASYMPTOTIC SOLUTIONS OF LINEAR ORDINARY 
DIFFERENTIAL EQUATIONS* 


By ROBERT L. EVANS (University of Minnesota) 


When asymptotic series are used in approximating the solutions of linear ordinary 
differential equations it is desirable to have an estimate of the accuracy of that approxi- 
mation. Older estimates of this kind are either difficult to obtain or of limited appli- 
cability. A simple method of making such estimates is presented here. It applies to linear 
second order differential equations, and for higher order differential equations it can 
be replaced by another and more laborious method.’ The results of this article can also 
be used to determine the number of terms of an asymptotic expansion giving the best 
approximation to the desired solution. 

The approximate solution of a given linear differential equation contains a fixed 
number of leading terms in an asymptotic expansion, and it differs from the correct 
solution by an error which, in its turn, is the solution of a second differential equation. 
This second differential equation is henceforth called the “related equation’ because it 
is intimately related to the given differential equation. Its appropriate particular solution 
is the error inherent in the approximate solution of the given equation.” However, the 
order of the related equation is one less than that of the given equation, so that when- 
ever the given equation is of second order the related equation is of first order. In this 
case of a first order linear related equation its appropriate particular solution (the error 
to be evaluated) can be obtained by merely multiplying by an integrating factor and 
then integrating. For the resulting error integral, an upper bound is easily estimated, 
but the obtained result does not in general agree with the widely-held belief that the 
size of the error is comparable with that of the first neglected term in the asymptotic 
expansion. In fact, it generally turns out to be a sum of multiples (or fractions) of more 
than one of the neglected terms. 

The description to follow in this article gives a general formulation of the error 
estimation, with a step-by-step parallel illustration of it in the treatment of a specific 
given differential equation. Equations of the general formulation are given numbers, 
and the corresponding equations of the illustration are given like numbers, together 
with the letter “a’’ after them—so, for example, (3) and (8a) would be corresponding 
general and specific equations, respectively. Also, the variables used in the specific 
example have a terminal subscript “a”. The selected illustrative example is shown later 
to have hypergeometric and Bessel subcases, but it is more general than either of these 
in order that it may also represent those cases in which the bound on the error contains 
more than one term. Finally, one of its subcases is used to provide an example in which 


*Received July 12, 1953; revised manuscript received Dec. 21, 1953. Presented in part to the Amer- 
ican Mathematical Society at Kingston, Ontario, Sept. 4, 1953. Developed under the Office of Ordnance 
tesearch contract DA-11-022-ord-489. 


‘Robert L. Evans, Proc. Am. Math. Soc., in press. 
2This scheme was conceived after reading Thorkell Thorkelsson’s work on “‘Divergent power series’’, 


Societas Scientiarum Islandica, 17, Reykjavik, 1934. 
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the smallest term of the asymptotic expansion is not the first term omitted by or in- 
cluded in, the best approximation of the desired solution. 
To specify the method of treatment, let 
K(w) = 0 (1) 
be the general given differential equation, in operator notation, and let one of its n 


linearly independent solutions be 
W(x) = ¢ y(x), (2) 
where y satisfies the differential equation” (the “given equation” of the preceding dis- 
cussion): 
Kle*y(x)] = L(y) = 2. [Ip.(z) ly” ” = 0, (3) 
— 4 


in which 


p(x) = Zz. b..* [vy = 0, 1, --- , n and each integer n(v) is finite] 


For a definite example consider 
L(y) = a*yh’ + (Aya + A,)2’y, + (Aga? + Aye + Ay. = 0 (A, #0). (3a) 


Assume that y(x) is such that 
y(z) ~ > (ele) |x~* (4) 


is a valid asymptotic representation in a particular sector S of | x | > x» . Analogously, 
assume that 


Y(t) ~ z=: [ca(a) Ja**~ “ ta) 
) (4a 


is the corresponding asymptotic expansion for the y,(x) of (3a). The powers p of (4) 
and p, of (4a) are later defined by (8) and (8a), respectively, and it is assumed that no 


logarithms need be introduced in (4). 
Before proceeding to further detail, let us outline the treatment by defining the 


approximate solutions of (3) and (3a) by 


u(x) = a [c(a) |x” , (5) 


a 


and 
1 


u(x) = z [ca(a) Ja** * 5a) 


respectively, where the fixed index N is arbitrary. With these quantities the associated 
errors of approximation are 


? p a . 
x)= y(@) — Ur) ~ 5 Cla) |x 6) 


and 
} . | { Ve z 4 
(t) = YAL) — U\t) ~ >. [Cl a) ja ; ba) 
a=N 
’The transformation from (1) to (3) was described in detail by C. E. Fabry (Thése, Paris, 1885) and 


in outline by the author in Ref. [1]. 
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respectively. In order to determine what differential equations have v(x) and v,(x) as 
their solutions we must next examine the nature of the c(a)’s and c,(a)’s. 

The c’s of (4) satisfy a recurrence relation—or a “serial relation’ in the terminology 
of Thorkelsson.* To describe the c(a)’s let 


m, = min [vy + nv)] = », + nr), 


Mz = max [vy + m(v)] = ve + mv), 


and 
m= ™M.— Mm, 
so that the recurrence or serial relation for the c(a)’s can be written in the form 
0 bh (p -atmy(p—-atm —1):-:(p—atm —n+t Dilc@ — m)} 


bh (p —-atm, — 1I(p—atm, — 2)---(p-—atm, —n) 


tb (p -atm — l\(p —a + m, — 2) -:> (p —-atm—n + 1) 
[ela — m, + 1] 


db, nop — a — Mm +m) +++ (p—a-—- m+tm —nt+rn + 1) 


+ other terms if v. is not unique 
‘[cla + m, — m,)], 


or, in reverse order and with simplifying notation, 


m 


b ds q,(a,p): [cla + m — v)] = 0, (7) 


v=0 


which is to be associated with the conditions determining the expansion in (4), namely 
that c(0 an arbitrary constant, and that c(— 1) = c(— 2) = --- = c(1 — m) = O. 
The analogous equation and conditions for the c,(@)’s are 


0 [4,(p, — a — 2) + Asg)[e(a + 2)] + As[c.(@)] 
+ [(p, — a — 1)(p, — a — 2) + Alp, — a — 1) + Ag Jle.(a + I] (7a) 


and 


If A, ~ 0, ¢.(0) is an arbitrary constant and c,(— 1) = 0, while 
if A, = 0, c,(0) is an arbitrary constant and (a + 1) in (7a) is replaced by a’. 


At this point one can define p as a root of the equation 
qo(—m, p) = 0 (8) 


and p, by 


p, = —A;/A,. (8a) 


‘Thorkel] Thorkelsson, “Serial relations and symbolic calculus”, Societas Scientiarum Islandica, 29, 


teykjavik (1951). 
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By using the results in the preceding paragraph we can next obtain the desired 
differential equations which are satisfied by v(x) and »v,(x). If we let 
| 


Cia) = ca +N), Ca) = ca + N), 
(9) 


o=ptN, and co.=ptN, 


then (6) and (6a) become 


v(xz) ~ > [C(@) x" 


and 
v(x) ~ p [(C.(a)]x"*~*, 


respectively. The insertion of substitutions (9) in (7) gives the recurrence relation that 

is satisfied by the C(a)’s and which corresponds to an nth order differential equation 
M(v) = 0, (10) 

where M(_) is a differential operator. The corresponding result for v,(x) is, from (9) 

and (7a), 

0 = [A,(o, — a — 2) + As + A,N][C.(@ + 2)] + A;[C.(@)] 


(¢, ~-a — l)(o, —a — 2) + Alo, —-a —1)+ Ay 
+ [C(a + 1)], 


+ 2N(o, —a — 1) + N(N — 1) + AN 
which corresponds to the differential equation 
0 = Miv,) = x’vs! + (Ayr + Ad)x’v, + (Aga? + Aye + As)0, 
+ 2N2°v. + No(A\x +N —1+ AQ, . (10a) 


In the final step of getting the desired differential equation for v(x) one combines (6) 


and (3) in the form 


L(y) = Lu + v) = Lu) + Lv) = 0 


to find, in view of (10), that 


M(v) — Liv) = 0 — [Lty) — LW) = Lt). (11) 


This linear differential equation in v(x) has an (n — 1)th order differential operator on 
the left-hand side and a known function of the last m terms of (5) (or all N terms if 
m > N) on the right-hand side. The discrepancy between u(x) and y(z) is that v(x) 
which is a particular solution of (11) and such that 

lim | 2*~'v(z) | = 0, 


2 @ 


in 8 
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in conformity with the representations in (4), (5) and (6). These results are more clearly 
illustrated by the specific example, for which (3a), (6a) and (10a) give 
M(v,) — Liv.) = 2N2°vo, + N2(A,x + N — 1+ A), 

= Lu) = (kz + A;s)Ty-1 + AsT y-2 , (11a) 
where 

k, = (p, + 1 — N\(p, — N) + Al — N), [in view of (8a)], 
Ty [c.(N — 1)]2°**'~*, 


and 


Tx-2 = [c(N — 2)]a’**?~*. 


Since (lla) is a first order linear differential equation for v,(x) the desired particular 


solution is 


v(x) = 2°~*~4”” exp (—A, 2/2) -I(2)/2N (12a) 


I(x) = | p(As+2Pe-N-3)/2 ayy (4 ¢/2)[(ki¢ + As)[e(N — 1)] + Astle(N — 2)]] dé, 


= £ 
g 
J —@/Ay 


where the path of integration is such as to make | exp (A, £/2) | nondecreasing. With 
this path of integration the exponential in the integrand may be replaced by exp (A, 2/2) 
to obtain the final error estimate, or bound, 

V — p, — 1)(N — p,) — AN — 1) 


FE > me N(A, — 2p, - mn D . T | 


4+ | A;/Nx(A, + 2p, — N — 1) |-| Ty-s | 
+ | A;/Nzx(Ao + 20. —-N+1)|-|Ty-2 |. (18a) 


Since the 7’s are the last terms of u,(x) this bound usually differs from the commonly 
guessed value, | Ty_, |, so that latter guess is usually erroneous. 

Discussion. The illustrative example has two or more well-known differential equa- 
tions as special cases. When A, = 1, As = k, Ay = 1/4 — m’, and A, = A; = 0, Eq. 


(3a) becomes the confluent hypergeometric equation and when A, = —2, A, = —A; = 
1, A, = — n’, and A; = 0, Eq. (3a) has as a solution the product of e* and the modified 


Bessel function of the second kind.° Since A; = 0 in both these subcases the last two 
terms in (13a) are then zero. 

When A, + 0 the recurrence relation (7a) is of second order and (12a) shows that 
error »(x) then depends on the last two terms in approximate solution u,(x), rather than 
just the last one. However, the terms involving A; [in (13a)] are often of smaller magni- 


tude than the other terms. 


‘=. T. Whittaker and G. N. Watson, A course of modern analysis, Cambridge Univ. Press and Mac- 
millan Co., N. Y., American ed., 1944, pp. 337 and 373. 
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In the usual use of asymptotic approximations it is a common practice to take N 
such that 7'y_, is the smallest term in the asymptotic series. In many cases it is not clear 
which of two successive integers is to be taken as N, and this uncertainty is said to reflect 
an error comparable to | 7'y_, |. Our (12a) and (13a) show that this earlier estimate may 
often be scarcely better than a crude guess and should be replaced by the present result. 

A further complication can arise when term magnitude | Ty_, | = | [e.(N — 1] 
z°**'-* | passes through more than one relative minimum with respect to N. In such 
cases it may happen that the term minimum of smallest magnitude does not occur at 
the N value giving the smallest error, v,(x). For instance, this arises in our example when 





A, = z=1,A; = — p, = 0.3, A, = 4.2 and A, = A; = 0. In this case | 7; | and | 7; | 
are relative minima of | 7'y_, |, and | 7, | is about 20 times greater than | 7’; |. But the 
error involved in taking (NV — 1) = 1 is less than a third of that involved in taking 
(N — 1) = 5 and less than a sixth of | 7’, | itself. 

One concludes, therefore, that for most accurate asymptotic approximation one must 
minimize the error, v,(x), with respect to the number of terms in u,(x) — N—for each 


specific differential equation and for each chosen value of the independent variable, z. 

This can now be done for second order differential equations by using (12a) or (13a) and 

for higher order differential equations by using the more complicated result of Ref. [1]. 
APPENDIX 


As a brief statement of the result in Ref. [1], as it applies to (3), if we define 


56= max _ [(degree of q,(a, p) ina — degree of go(p, a) ina)/r], 
if 6 < 1 (which, if not already satisfied, can be brought about by substituting « = 2” 
with M a sufficiently large positive integer), and if | w | = 1, then for any N 2 2 
rs) 8 6 
xz v(x) = > Zz A,.a(w)” “[e(u + pi / I] (x + pw), 
B=N-1 p=N-1 u=0 
[where Ay, = A,.,-1 = 0, A:,, = 0, and, for B > 1, 


Aye = (8B — LNA, e-1) + Aww 1 B 1] 


in| z| > 2, and|argz — argw| S 2/2. Incidentally, arg w is arbitrary. 


WAVE PROPAGATION IN A VISCO-ELASTIC MEDIUM* 


By E. J. SCOTT (University of Illinois) 


1. Introduction. The use of mechanical models to describe the behavior of materials 
is well known [1]. By combining a Maxwell unit, consisting of a spring and dashpot 
connected in series, with a Voigt unit, consisting of a spring and dashpot connected in 
parallel, in various ways, it is possible to describe adequately in many cases the complex 
behavior of real materials. The stress-strain laws corresponding to a Voigt unit and a 
four parameter model consisting of a Maxwell unit combined in series with a Voigt 





*Received Sept. 3, 1953. 
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unit have been used to study the physical properties of high polymers and wave propa- 
gation in a visco-elastic medium [2]. 

It is the purpose of this article to consider the propagation of longitudinal waves in a 
medium confined between the two parallel planes x = 0 and x = I (Fig. 1) and whose 





Fria. 1 
physical behavior is representable by the linear mechanical model shown in Fig. 2. 


The particular case of two Maxwell units coupled in parallel will be considered in detail 
because this model seems capable of describing some aspects of creep. 




















2. Derivation of the partial differential equation. We shall first derive the stress- 
strain law that obtains for the linear model of Fig. 2. From elementary considerations, 


we have the following equations: 


K ;[e,(x, t) — €:,(x, t)] = Exe(a, 0, pw 5S. ore ee (1) 
me, (x, t) = o(2z, t) — Zz K,le(a, t) — «..(2, 0], (2) 


where x is a parameter and the subscripts involving ¢ denote differentiation with respect 


to time. By elimination, we obtain 


o(x, t) = me,,(x, t) + =. Exe(z, t). (3) 


i=1 
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Solving Eq. (1), subject to the condition e;(z, 0) = 0, for e,;(z, ¢) and substituting 
in (3), we find that the stress-strain law for the model under consideration is 


a(x, t) = me,,(x, t) + pe BE. | ex, rT) exp [—yu,(t — 7)] dr, (4) 


OG 
where pn; = E,/K; 
In the slab of thickness | (Fig. 1), let u(x, t) represent the displacement of the section 
x at time ¢ and let p be the density of the unstrained medium. Then 


o,(z, t) = pu,2, 0, (5) 


ea, t) = uz, o), (6) 


where the subscripts denote partial differentiation with respect to the corresponding 
variable. 

Eliminating o and e we obtain the following integro-differential equation for the 
propagation of longitudinal waves in the medium under consideration 


pu,,(2, t MUzri(2, t) + ¥ E; Usey\L, T) CXDP [—p(t— t) | dr. (7) 
1 l “0 


3. Solution. In what follows we shall assume that the material is initially unstrained 


and at rest, so that 


u(x, 0) = u,(z, 0) = O. O<2< J. (8) 
Applying the Laplace transform 
u*(z, p) = L, {u(z, )} = | e ‘u(x, t) dt 
to Eq. (7) and making use of conditions (8), we obtain 
d-u*(x, p) du*(x,p) < E; 
ppu*(x, p) = mp — a ——— a 
f f f i er 2, Pt ui’ 
from which 
d-u*(zx, p) opu*(x, p) 
= = th 8 = 0. (9) 
ax . , 
mp + > E; (p + Hi) 
Rad | 
The solution of this equation suitable for the finite region being considered is 
u*(x, p) = c, cosh Zi” + co sinh ry lal (10) 
where Z, = pp/[mp + Dui-1 E:/(p + us)]- 
Let us now suppose that the section x = | is fixed and the section x = 0 is given an 
impulsive velocity V, which is then maintained, i.e., let 
u(0,=Vo., ulj=0, (t>0), (11) 


whence 


u*(0, p) = Vo/p u*(l, p) = 0. (12) 





r 
> 
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These conditions when applied to Eq. (10) yield the equations 
V./p =a, 0 = ¢, cosh 1Z'” + ¢, sinh 1Z'”. (13) 
The values of c, and c, obtained from Eqs. (13), when they are substituted in Eq. 
(10) and the result is simplified, give 


V, sinh (1 — 2)Z,”” 











u*(z, p) = ——s wei : 1 
(x, P) p sinh IZ)” Sas 
whence, by the complex inversion theorem, 
Vo f°" _,. sinh (l — 2)Zi” 
At, GO = = | oo ad 15 
) Qi Jie © pp sinh IZ? (15) 


4. A special case. We shall consider the solution, in detail, for the special case of 
a model for which the number of Maxwell units coupled in parallel is two and m = 0. 
In this instance Eq. (15) becomes 
V, perio . Si l = — 
u(x, ) = = | p Ca — 20s ay (16) 


<a e€ on 
2m Je-ic p sinh IZ, 





where Z, = pp(p + w:)(p + m2)/[(E, + E2)p + Eis + Exum]. 
The integrand is a single-valued function of p with a double pole at p = 0 and poles 
at those values of p which satisfy the cubic equation 


2. S 
uw 2 


2. 3 
p' + op" + ("os + oi)p + ot =0, n=1,2,3,---, (17) 


where o; = 4; + uo, 02 = Mita, 03; = E, + E., 04 = E,u2 + Eu, . Since the coefficients 
of Eq. (17) are all positive, its roots, if they are all real, will all be negative. Further- 
more, since 
(ae 2) nn” , — 
O1\ sa 03 + 2 > “a Oe 5 n= 1,2,3,---, 
l"p / lp 
if Eq. (17) has one negative and two conjugate complex roots, the real parts of the 
complex roots are negative [3]. 

Making the substitution p = y — o,/3 in Eq. (17), we find that 


nr 


x 2 3 2/ 
+ | (ot — o,0;/3) +5501 — 2.01/3 | = 0, n= 1,2,3,°°-. (18) 
ad 


Case I. If n’x’e;/(l’p) + 03 > o:/3, n = 1, 2, 3, --- , then the coefficient of y in 
Eq. (18) is positive and Eq. (18) has one real root and two complex roots. Consequently 
the roots of Eq. (17) can be represented thus 


p= —a,, p= -b+t,, p=-—-b— t,, 
where a, > 0, b, > O and n = 1, 2, 3, --- 
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The residue of the integrand of Eq. (16) at the double pole p = 0 is 


— 2 ain ee ie De 
l - 2 t+ poox(1 _ IX - 21) 
l 6lo; 
and that at p = —a, is 
e”' sinh (1 — x)Z; | ("2 Qnr(o, — o3a,)e °" sin [nr(l — x)/T] 
er ome _ —_— _—_— } 5) a = _~ ~ = 77 — se eee © 
p d(sinh 1Z,")/dp J, a,l’[p(3a, — 20,4, + 02) + on /T| 


Furthermore, the residue at p = —b, + tc, 1s 


e”' sinh (1 — 2)Z; | 


p d(sinh 1Z,’*)/dp 
2nr[(A,C, + B,D,) + (AD, — B,C,)] sin [nw(l — x)/l] exp (—}, + 2,)t 


PA? + Bb — cy + 40] 


= (—] 


and that at p = —b, — ic, is 


(—1)’ 2nn[(A,C, + B,D,) — (A.D, — B,C,)] sin [nx(l — «)/l) exp (—b, — tc,)t 


[?( 42 + Be)((b2 — c2)’ + 4d5c;] ; 


where 


A, = 3p(b; — c,) — 2e,pb, + po. — o3n'r /T, 


C, = o1b5 — bia, — osc, — 0,C,03 , 
D. = Qb.C.0% — 636, — b*c,.03 : 
Therefore, 
( 
u(x, t) = V,\(l — x)t/l + pozx(l — x)(x — 21)/(6le;) 
a 2m Yi (—4)"*} nor — o3d,) sin [na(l — 2) Ue ae 
; = a,|p(3a%, — 20,a, + o2) + on /T] 
* Lr > (—p n(M; + N;) “¢ | bnt C08 (c t - $,) sin Ina(? — x) yy (19) 
[? : (A + B;)[(b, — ¢.)” + 46:c;, 
where MVM, = A,C, + B,D, , N, = A.D, — B,C, and tan ¢ = N,/M, . 
Case II. Suppose that n’x’o,/(l’p) + 02 < 01/3 for n = 1, 2, --- , k, where k is 
the largest integer for which the inequality holds. Then Eq. (18) may have three real 
roots and Eq. (17) three real negative roots. For n = k + 1, k + 2, --- , however, 


n*x’c,/(?p) + o2 > o;/3 and the solution will have the terms given in Eq. (19). There- 
fore, if the first inequality mentioned above obtains for n < k and the second for n > k, 
then the solution will contain the terms displayed in Eq. (19) plus (possibly) three 
summations of terms like those occurring in the first summation of Eq. (19). 

An examination of the solution reveals that ast —~ © the transient terms following 
the summation signs become negligible and the oscillations tend to a state characterized 
by the first two terms of Eq. (19). 
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5. Discussion of a numerical example. Figures 3 and 4 show the plots of u(z, t)/Vo 
and u,(x, t)/V> versus t, respectively, for different cross-sectional values z of a medium 
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characterized by two Maxwell units coupled in parallel and such that m = 0, EZ, = E, = 
K, = K, = p = 1 = 1. Both graphs indicate the damped oscillatory character of the 
displacement as well as the velocity of the cross-sections between x = 0 and x = 1 
of the medium. As t — © the oscillations die out and the displacement of a section 
increases uniformly while its velocity becomes constant. While the graphs depict the 
general character of the motion, computations were carried out for a limited number 
of points and exact values are not implied at all points. 
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NOTE ON TAYLOR INSTABILITY* 
By GARRETT BIRKHOFF (Harvard University) 


1. Qualitative discussion. In a well-known paper [3], Sir Geoffrey Taylor has 
discussed the stability under normal acceleration of a plane interface separating two 
fluids of different density. His main conclusion (reached several years earlier) is now 
classic: the interface is unstable when the light fluid is accelerated towards the dense fluid 
and (presumably) stable when the reverse holds. This conclusion has important applica- 
tions to gas-filled underwater explosion bubbles. 

Its applicability to small vapor-filled cavities is however less clear. The stabilizing 
role of surface tension is known’ to be important, and Binnie [1] has suggested that 
surface tension may even be sufficient to compensate for Taylor instability. 

The purpose of this note is to show that, in spite of the fact that the denser liquid is 
being accelerated towards the lighter vapor, collapsing bubbles are unstable, and that this 
result is unaffected by surface tension (though it may be affected by viscosity or thermo- 
dynamic considerations). The proof of this fact depends on a consideration of the 
stability of differential equations near regular singular points: the instability is algebraic, 
and not of the exponential type usually considered. 

2. Negative damping. The formulas underlying the perturbation theory of collap- 
sing spherical cavities are easily found’. Let the cavity radius b(t) be given as a function 
of time, and the interface expressed in spherical coordinates in terms of Legendre poly- 


nomials by 
r = b(t) + >> b,()P,(cos ¢). (1) 
h=1 


Supposing the b,(¢) small, and neglecting gravity and surface tension, the condition of 


constant internal pressure givesT 
bb,’ + 3b'b, — (h — 1)b'b, = 0. (2) 


The same formula applies to any surface harmonic of order h. 


*Received Oct. 19, 1953. 
IThis was observed independently in 1951 by R. H. Pennington and R. Bellman at Princeton, and 
by R. L. Ingraham and the author, but not published. 

2Formulas (1)-(2) were derived by W. G. Penney and A. T. Price, British Report SW-27 (1942); see 
R. H. Cole, Underwater explosions, Princeton, 1948, p. 311. See also [1, Formulas (2)-(3)]. 

tFor typographical reasons, the dots indicating time derivatives are set as superscripts. 
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The facts about stability are suggested by the ordinary stability test. Near b = 0, 
b’° < 0, and so the coefficient of b, in (2) is positive—i.e., we have ““Taylor stability” in 
the usual sense*®. However, the coefficient of b; is negative in the collapse phase, and so 
we have “negative damping”. This suggests that vapor-filled bubbles are unstable 





during collapse. 

However, it is more satisfactory to get quantitative results; this we will now do. As 
a by-product, it will appear that the ordinary stability test is not conclusive (see Theorem 
2). 

3. Asymptotic analysis. Since the kinetic energy is proportional to b*(b’)’, [2, p. 122, 
Eq. (6)], which approaches a constant value as b | 0, evidently b'*°db and dé are 
asymptotically proportional to each other, whence ¢ « b*’ and b « ¢°**. Substituting 
in (2), we get the asymptotic equation 


7 Sy ee 24(h — 1) 
b, + i : bi, + a, e 


bh, = 0, near b = 0. (3) 
The sign of the coefficient of b; depends on whether the expansion or collapse phase is 
involved (i.e., on whether ¢ > 0,b° > Oorb’ < 0,t < 0). 

Equation (3) has a regular singular point* at ¢ = 0; the indicial equation, corre- 
sponding to possible exponents in asymptotic complex solutions 


by, = (Co +t + eet” + ---) (4) 
of (3), is 
ala — 1) + 1.2a+4+ .24(h — 1) = 0. 
This has the conjugate complex roots 
a=—1+ V0l — 24h — 1), (5) 


which corresponds to b, « t°' « b7'“*. We conclude 


THEOREM 1. The amplitude of any perturbation is inversely proportional to the 
fourth root of the radius, as the radius of a collapsing spherical cavity shrinks to zero. 

In particular, we have stability during expansion, and instability during collapse. 

A consideration of the imaginary part of a shows that the phase of the oscillation, 
which is approximately the same as that of 


undergoes a half-oscillation when b changes by a factor of #/1.25 V h — 122.5/Vh — 1. 

Of course, when b,/b exceeds 0.2 or so, linear perturbation theory is no longer appli- 
cable, and Theorem 1 ceases to be physically applicable. (The same limitation applies 
to Binnie’s discussion. ) 

4. Physical generalization. The preceding analysis can be extended by a few re- 
marks to other physical situations. 

This seems to be the only sense considered in [1]; see [1, formula (14)]. 

‘For regular singular points and their indicia] equations, see L. R. Ford, Differential equations, p. 
184 ff. The inadequacy of Binnie’s approximations should now be evident. 








308 NOTES [Vol. XII, No. 3 

Remark 1. Surface tension has a negligible effect asymptotically. Thus its asymptotic 
contribution to the available energy is negligible, so that it does not affect the formula 
ba« #* b « bb” &« b-*. Now, writing (2) in the alternative form 


bi’ + 3Cb-?*b; + 1.5C(h — 1d, = 0, (3’) 


we see that (i) surface tension does not affect the coefficient of b; , and (ii) since surface 
tension alone causes oscillations of a period proportional to b’**, the contribution of 
0(b-*) to the coefficient of b, is negligible. 

Independently, a simple consideration of reversibility (the invariance of (2) under 
b —— b, i — — #) shows that positive stability in the expansion phase must be matched 
by instability during collapse. 

Remark 2. Similar calculations hold for collapsing cylindrical cavities; the 3 in (2) 
must be changed to 2, and the p,(cos ¢) in (1) to cos hg. The kinetic energy varies ac- 
cording to a less simple law’, namely b’(b’)’ In b. Neglecting the logarithmic term for 
simplicity, we set b’ < b™’, orb « ¢°’*. This gives 


(h — 1) 
4t° 





bv + 5 bi + b, = 0. (6) 
The indicial equation is a* + (h — 1)/4 = 0. Hence, to this approximation, the amplitude 
of perturbations neither grows nor decreases. In this sense, we have neutral stability to 
a first approximation; in the sense of relative amplitude 6,/b, we have instability during 
collapse, as in the spherical case. Finally, we have a half-oscillation in b, when b changes 
by a factor of 2r/Vh — 1. 

5. General theorem. More generally, one can analyse the stability of solutions of 
ordinary linear differential equations near regular singular points, by writing down the 
dominant terms. Thus consider 





. d’x — d”" ‘x as 
t ‘Ta a,t FT aa +---+a,z = 0. (7) 


The following result is self-evident, in view of (4). 

THEOREM 2. Stability occurs as ¢ 7 0 if the real parts of all roots of the indicial 
equation‘ of (7) are positive; instability occurs if any of them is negative. 

Applying the usual Routh-Hurwitz criteria’ for this to the indicial equation, we 
derive as corollaries the following special tests for stability of (7). 

Case n = 1. One condition, a, < 0. 

Case n = 2. Two conditions, a, < 1 and a, > 0. 

Case n = 3. Three conditions: a, < 3, a; < 0, (3 —a,)(2 — a, +a.) > —a;. 

Case n = 4. Four conditions: a, < 6, a, > 0, a, > 3a, — 11, and 


(6 — a,)(az — 3a, + 11)(6 — 2a, + a, — a;) > (6 — 2a, + a, — a3)? + (6 — a,)’a,. 


Remark. Note that, even when n = 2, the tests on the coefficients of (7) are more 
subtle than the usual tests. 


‘T. E. Sterne, J. Appl. Phys. 21, 73-74 (1950). 
‘See J. V. Uspensky, Theory of equations, New York, 1949, Appendix III; by Theorem 2, — a must 


replace a in the usual test. 
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Coro.iary. Suppose that 
bb, + Ab'b, + Cb''b, = 0, (8) 


and that b « |¢|°asb +0,t—0. Then: 
(i) If C8(@ — 1) < 0, almost all b,(¢) increase without limit ast T 0 (ort | 0). 
(ii) If C8(8 — 1) > O and AB > 1, then all b,(¢) increase without limit ast 7 0 
(ord | 0). 
(iii) If C8(8 — 1) > O and AB < 1, then all b,(¢) tend to zeroast T O(ort | 0). 
Proof. Substituting b = Ke?’ in (8) and simplifying, we get 





tb,’ + ABth, + CA(B — 1)b, = 0. (9) 
This has a regular singular point at ¢ = 0, with the indicial equation 
a’ + (AB — l)a + CBB — 1) = 0. (10) 
The two roots a, , a, of this are 
1 — AB 





Qa; = 


1 3 
—~ ws V/(AB — 1)? — 4CB(6 — 1). 


In Case (i), a, and a, have opposite sign, and so the general solution b, = c,t** + c,t** 
of (8) (logarithmic terms* do not alter this) increases without limit as ¢ f 0. In Cases 
(ii)-(iii), a, and a, both have the same sign as 1 — A§, whence the conclusions are still 


obvious. 
I wish to acknowledge my indebtedness to Dr. R. L. Ingraham for helpful suggestions 


about this paper. 
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ELECTROMAGNETIC WAVE PROPAGATION IN BOUNDED ELECTRON BEAMS* 
By PHILIP PARZEN (New York University) 


The linearized Maxwell’s equations in MKS units for the field quantities in electron 
beams are: 


V XE = —jonH, (1) 

V XH =J + cE + jwcE, (2) 

V-J + jwop = 0, (3) 

joV + Vo VV +V-VV. = — [E+V XH, + V, X Hl, (4) 
J _ poV + pVo . (5) 


Received October 20, 1953. 

*This work was performed at Washington Square College of Arts and Science, New York University, 
and was supported in part by Contract AF 19(122)-1, with the U.S. Air Force, through sponsorship 
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Here all field quantities are the sum of a time independent and a sinusoidally time 
varying term. E, H, J, V, » denote the time varying electric field, magnetic field, electron 
current density, velocity and charge density, respectively. The subscript “‘0’’ denotes 
the corresponding time-independent quantities. €) , uo , ¢, — e and m are the dielectric 
constant, permeability, conductivity of the medium, charge and mass of the electron, 
respectively. The problem is, given the time independent quantities, to solve Eqs. 
(1)-(5) for the time varying quantities subject to boundary conditions chosen in accord- 
ance with the following uniqueness theorem to be proven. 

In the case that J is parallel to V, the time varying field quantities within a surface 
S are uniquely determined by specifying on the boundary either E X n or H X n and 
V,:V or J-n. n is the unit normal to the surface. The proof is as follows: From Eqs. 
(1) and (2), one gets, denoting conjugates by an asterisk, 


[xz <x H*-dS = al [ (upH-H* — «,E-E*) ar | + / oE-E* dr + [ J*-Ear. (6) 


| 
* 


Using the vector equation identity, 
V(vo-V) =V-VV+VVV+V-VXVt+VVXYV~. 


Equation (4) gives: 


joV + V(V,-V) = —(E+VxH+VxXHI-VWxXVxXV-VxVxXxV.,. 


t 
in 
Since J and V and hence J, and V, are all parallel to each other, 


” 


JE = 7 [joV-J* + J*-VV-V)]. 


Using Eqs. (3) and (5), the following conservation law is obtained: 


“ 


; F m i 
iol | uoH:-H* — ¢E-E* + pv-V* | dr + | oE-E* dr 
. ft 
| 
: m -_ 
+ | | <x H* + wv,-W)y* |-as 0. (7) 
« t 
The proof of the uniqueness theorem follows immediately by taking the real part of 
Eq. (7). This conservation equation can also be used to define impedances and scattering 
matrices [1]. The condition for oscillations in a resonant cavity is, 


| E H-H* — «E-E*+™ pv-V" | dr = 0. 8) 


The case J parallel to V occurs, for example, when H, ©. In this case J, V and Hy, are 
parallel to each other. 

A solution will now be given for the modes of propagation in a right cylindrical 
waveguide (axis that of the Z-axis) which is completely filled with an electron beam 
which is constrained to move solely along the Z axis (HW) = @ and along Z-axis) and 
whose time independent velocity is an arbitrary function of Z. According to the unique- 
ness theorem, the fields in a volume G bounded by the planes Z = Z, , Z, , and the wave- 
guide surface is uniquely determined by specifying at Z = Z,; E X Zand J = E,(0) 
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and J(0), respectively; and similarly at Z = Z, . As in the case without an electron 
beam [1], the field can be expanded in terms of transverse magnetic (7'M) and transverse 
electric (TE) modes whose coefficients are functions of Z. The transverse dependence 
of these modes are exactly the same as in the space charge-free case. The longitudinal 
dependence of the 7E modes are also the same as in the space charge-free case. The 
longitudinal dependence of the coefficients of the 7M modes is given as follows. Let 
J,V,V.,J., Hand H denote the longitudinal component of these quantities, and the 
subscript Z denote differentiation with respect to Z. Then, for any mode, 


Vidaz + [2jaVa + 3ViVoz]I 2 + [2joVoVoz — w' Vol = = jok, (9) 
Err + (k° — TE = oi. (kJ + J22), (10) 

-TE XZ = = is ~ th, (11) 

HX Z = — - ((E X Z), — E]. (12) 


T, is the transverse propagation number of the mode and k? = w’yo¢ . The coefficients 
of the 7M modes are uniquely determined by the boundary conditions previously 
specified which is equivalent to specifying J and Jz + jwe/z at Z = Z, , Z, . Equations 
(9) and (10) constitute a two point boundary condition linear differential system [2] 
which possesses a unique solution if the system with null boundary conditions has only 
the trivial solution J = E = 0. This will be the case according to the uniqueness theorem 
which implies that zero boundary excitation leads to zero fields in the interior. It is also 
easy to prove that the fields are uniquely determined by specifying J, EZ, Jz and Ez 
or J, Jz,E X nandH X nat Z = Z, and nothing need be specified at Z = Z, . 

To obtain explicit solutions, Eqs. (9) and (10) are to be solved for J and E subject 
to initial values of J, E, Jz and E;z at Z = Z, . Writing the equations in matrix form, 











Uz, =AU; (13) 
| J 
Jz 
U = column matrix - U = U(Z,) at Z=2,; (14) 
jwoE 
Ljwk z 
0 1 0 0 
Yo + J2Y02z —J2y7o + 3¥02z/y0 y'¥o 0 
matrix A = > (15) 
0 0 0 1 
—k? — ¥5 — j2voz j2v0 — By¥02z/¥o —k’ +Ti- yy 0 
w® 2 é . 
Tee 5 7; 2° 
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As suggested by Ramo [3] in the case of Voz = 0, a solution is sought which may be 
expressed as the sum of four exponential functions whose coefficients are uniformly 
convergent power series in terms of y* and Voz . A convenient method to accomplish 
this is the matrizant scheme [2, 4]. The solution is, 


Z 
U(Z) = 2 (ASU(Z,). (16) 
Zz , 
Q {A} = matrizant of A 
ia aZ ~Z Xs 
ade ihe 2 | A(X,) dX, + | A(X.) dX, | A(X,) dX, + -+:. 
“Ze “Ze “Ze 


This solution is converted to exponential form by transforming A or a portion of A into 
diagonal form [5]. Only the final result will be given here [6]. 


U = PiD3[Py"(Zo) + WPr'(Zo) + 6(Zo) + ¥o(Zo)]U(Z.). (17) 


D is a diagonal matrix = diagonal (d;), 7 = 1, 2, 3, 4. 


aZ 
‘D? = diagonal (exp d; ax). 
“Ze 
aZ 
U. = > b (Py'(Zo) + WPr"(Zo) + 6(Zo) + ¥O(Zo)) ix[Ui(Zo) Pri exp | d, ax | 
i k Ze 
(18) 
The coefficients of the exponentials are uniformly convergent series for small values of y’. 
Conclusions. A uniqueness theorem for the sinusoidal fields in electron beams has 
been proved for the case where the current density is parallel to the electron velocity. 
The modes of propagation in a right cylindrical waveguide completely filled with an 
electron beam with an arbitrary longitudinal time independent velocity has been worked 


out. 
Acknowledgement. The author wishes to thank Professor Friedman for his aid. 
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A NOTE ON THE NUMERICAL SOLUTION OF THE EQUATION 
x” + ax” + b = 0* 
By TING KUAN SHU-CHUANG (United Nations) 


In this equation, a and b are real and not equal to zero; m and n are positive integers, 
m>n.™ 

We give a method of transforming this equation which permits us to obtain the 
complex non-real roots by solving independently for the argument 6 and the modulus r. 

Substituting x = r exp (7@) in the above equation yields, for the real and imaginary 


parts, 


r™ cos m@ + ar" cosn@ + b= 0 (1) 
and 
r™ sin m@ + ar’ sinné = 0. (2) 
Solving (2) for r, 
1 
asin n6@\" " 
p= (-s8nnty” oo. (3) 
sin mé 
Substituting (3) in (1), we have 
asin n0\""" asin n@\”" " 
(- : ane) cos mé + a( - —— cosn@ = —b. 
sin m@ sin m@ 
Factoring, 
sin n6é \"~" sin n@ rs 
(— _— ) (cos n@ — -—— cos mo) = —b/a”™™. 
sin mé sin m@ 
Simplifying, we conclude 
sin n@ \"~"[ sin (m — n)@ an 
(—sinae) si (m — s)e = -b/a. (4) 
sin mé sin m@ 


We have therefore an equation in @ alone. Clearly it can be solved by repeated 
substitution, and then r can be obtained from (38). 
The fact that r is positive can be used conveniently to limit the range which we must 
test for 6. Equation (3) implies 
yn-n _ _ a SIN NO (5) 
sin mé 
Thus sin n@ and sin mé@ have opposite signs if and only if a is positive. On the other hand, 
if we multiply (4) by a’’"~", it becomes 
{rnin = 98] _ 2. « 
sin mé a 


*Received October 28, 1953. 
**The restriction on m and n is not necessary but the use of more general m and n requires a defi- 
nition of the corresponding powers of a complex number. This would complicate the discussion. 
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Hence sin (m — n)@ and sin mé@ have the same sign if and only if a and b have opposite 
signs. These conditions reduce considerably the amount of test work for @, as will be 
shown in the example. 


6 


Example. x + 2 +1 =0. 
According to (3) and (4), we write 


sin 66\'”? 
a = € ) (7) 


sin 86/ 
and 
( sin 66) bes ae) 
~ sin 86/ \sin 86/ =a, 
whence 
sin® 6@ sin 20 = sin‘ 88. 


The limiting condition for the range of @ is that both sin 26 and sin 60 should be of 
opposite signs from sin 86. Thus we need only to consider the following ranges for @. 


(a) 22.5° — 30.0° (e) 202.5° — 210.0° 
(b) 67.5° — 90.0° (f) 247.5° — 270.0° 
(c) 90.0° — 112.5° (g) 270.0° — 292.5° 
(d) 150.0° — 157.5° (h) 330.0° — 337.5° 


These ranges can be paired off as (a), (e); (b), (f); (c), (g); and (d), (h); or, as (a), 
(d); (b), (ce); (e), (h) and (f), (g). For the first group of pairs, the values of the sines of 
26, 6@ and 8@ are identical for each pair, respectively, whereas for the second group, 
they are symmetrical and have opposite signs for each pair, respectively. Thus once the 
test for @ is completed for any two ranges not belonging to the same pair, the other 
values of @ follow automatically. 

Suppose we take (a) and (b). Owing to the periodicity of the sine function, it can be 
determined by inspection that @ lies within 25.5° — 30.0° and 72.0° — 77.5°. With the 
aid of a sine table, the respective ranges can be narrowed down further by inspection to 
25.5° — 26.5° and 74.0° — 75.0°. From this point the location of @ can be carried on by 


mt 


computation. 
The approximate values of @ are: 
26.05077° 206.05077° 
74.70037° 254.70037° 
105.29963° 285.29963° 
153.94923° 333.94923° 


Substituting in (7), we obtain the respective approximate values for r: 


91911 91911 
1.08800 1.08800 
1.08800 1.08800 

91911 91911 





e 
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Therefore the respective approximate values for x are: 


£82573 + 403647 — 82573 — .403647 
.28709 + 1.049447 — .28709 — 1.049447 
— .28709 + 1.049447 .28709 — 1.049447 
82573 + 403647 82573 — 403647 


NOTE ON THE CIRCLE THEOREM OF HYDRODYNAMICS* 
By I. LEVIN (National Bureau of Standards) 


The circle theorem [1, 2] is concerned with the irrotational, two dimensional flow of 
an incompressible, inviscid fluid in the z plane. Let f(z) be the complex potential of flow. 
Then if there are no singularities within a radius, a, of the origin and no rigid boundaries 
in the plane, the appropriate flow function after introducing a circle | z | = a about the 
origin is given by g(z) = f(z) + f*(a°/z) (where * denotes conjugation). 

The purpose of this note is to show that the restriction on rigid boundaries may be 
somewhat relaxed. 

That some restriction is necessary may be seen from the following example. Consider 


the uniform flow g(w) = w. From the circle theorem, h(w) = w + 1/w represents the 
flow with a unit circle about the origin. By the transformation z = w + 31, we see that 
H(z > — 3i + (z — 37)~* represents the flow about the circle centered at 37. Then if a 


second unit cirele is introduced about the origin, F(z) = 2 + (z — 31) +2° 4+ (27 ote 
3i)~' should give the flow for the two circles. It can easily be verified that while z = e*° 
is a streamline, z e* + 32 is not. 

A rigid boundary will be called generally admissible if for every flow f(z) having 
this boundary as a streamline and satisfying the remaining conditions of the circle 
theorem, the function g(z) = f(z) + f*(1/z) has both the unit circle and the boundary 
as streamlines. The boundary will be called conditionally admissible if this is true only 
for certain complex potentials f(z). 

To determine a necessary and sufficient condition for rigid boundaries to be generaily 
admissible, consider any complex potential f(z) = a(x, y) + 78(x, y). Then 


[f(z)]* = a(x, y) — iB(a, y) = f*(e*) = f*(x — wy), 


x Yy x y x y 


ve g) == f* a 7 5} = ae a! =] — Bla 3) 2 2), 
rae J (2 + ¥" z+ 3) ae +y z« + ) ia +y x + -) 


, _ j x y 

S{g(z)] = dl f(z) f*(1/z)] = B(x, y) — ( — s ¥,). 

i[g(z)] [fe +. | = A(x, y) — Ala ge a ae 

Clearly if x7 + y° = 1, s[g(z)] = 0, and hence the boundary of the unit circle is a stream- 
line. Suppose C: x = 2(r), y = y(t), To < 7 < 71, is a rigid boundary. Then 9[f(z)] = 
Bla(r), y(r)] = constant for tr, < +r < 7, . In order that this rigid boundary be admissible, 


7 


it is necessary and sufficient that 9[g(z)] be constant along C. Hence B[x/(a* + y’), 


*Received November 27, 1953. The preparation of this paper was sponsored (in part) by the Office 
of Naval Research, USN. 
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y/(x’ + y*)] = constant along C. That is, the inversion of the boundary with respect to 
the circle must also be a streamline of the original flow f(z). To be generally admissible 
this inversion of the boundary must be a streamline for every flow hence must itself be 
a rigid boundary. Thus a necessary and sufficient condition that rigid boundaries be 
generally admissible is that they be mapped into themselves under inversion with respect 
to the circle. A given boundary can be made conditionally admissible by obtaining the 
flow for this boundary together with its inversion and then regarding the inversion as 
only a streamline and not a boundary. 

In particular, any plane barrier along y = cx, a < x < b, is generally admissible 
subject only to the condition that ab = (1 + c’)~*. As an example of the application of 
the circle theorem under this modified restriction on rigid boundaries, consider 

9 8. | 1\ ]}}2 
f(@@) = 5 cos 6[42 — 5] — i ssin 6 (z — 2)(c _ 1) ] 
3 é 2 
representing flow at an angle @ about the plane boundary 1/2 < x < 2. Then under the 
transformation g(z) = f(z) + f*(1/z) it can be shown that both the unit circle and the 
plane boundary are streamlines. Any portion of this plane boundary is conditionally 
admissible with respect to this complex potential f(z). 

Acknowledgment. The author should like to express his gratitude to Professor R. 

Redheffer for helpful suggestions and criticism. 
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A NOTE ON A PAPER BY G. C. McVITTIE* 
By G. B. WHITHAM (The University of Manchester) 


In a recent paper** in the Quarterly of Applied Mathematics G. C. McVittie uses 
Einstein’s equations in general relativity to derive certain solutions of the classical 
equations of continuity and momentum for the compressible flow of a fluid when heat 
conduction and viscosity are neglected. Explicit expressions (which always satisfy 
these equations) are obtained for the density p, pressure p and velocity components U;,; 
in terms of arbitrary functions ¢, ¢; , ¢2 , ¢; of the time 7 and space coordinates X, ; 
the ¢’s are not independent however and must satisfy certain (differential) consistency 
equations. The equations of continuity and momentum are (using the double suffix 


summation convention) 








9p 5 9 ry _ 
oT ax, (pl ) a7 0, (1) 
au, : ou) — “ 
of oT a ihed, OX; + \ 0; (2) 


*Received December 2, 1953. 
**G. C. McVittie, “A method of solution of the equations of classical gas dynamics, using Einstein’s 


equations’’, Quart. of Appl. Math., 11, No. 3 (1953). 











1954] G. B. WHITHAM 317 


MeVittie’s solutions are 





= —V%, (4) 
pe $43 1 E(ou-3 oe) +5 or $7) /v%. “ 


i=-1 





where, if l, m, n is a cyclic a of 1, 2, 3, 





























a’¢: ‘/ v6 (6) 
OX, OXn = ox es 7 oT : 
and 
d°b> , obs ( ve) 1 _ Hb , ibs ( ve) eo 
ax: 7 ax: + aX, dT/ Vo aX; + oxi + ax, dT/ Ve 
_ 9) , de ( ve) eS 
~ @X? + ax? + aX, oT/ Vo’ @) 


In any problem ¢, , ¢2 , ¢; must be eliminated from (6) and (7) so that p, p, U, are ex- 
pressed in terms of ¢ alone; usually ¢ would then have to satisfy some energy relation. 
The purpose of this note is to point out two results which would be of considerable 
value to anyone using this solution: (i) the elimination of ¢, , ¢2 , ¢; can be carried out 
explicitly, and (ii) the results can be obtained in a few lines (and in fact more generally) 
directly from (1) and (2). This latter derivation is set out first. 
It is clear that (1) can be satisfied by the introduction of an arbitrary vector A, such 


that 


Q 





— 
piigues ax,’ pu, i ee aT ° (8) 


Equation (2) then expresses dp/@X, in terms of A; . The simplest form is obtained if 
(1) is multiplied by U, and added to (2) to give 














ee 
ax, = (eU J+ 3 (pU.U;), (9) 
whence, on substitution from (8), 
|) ee OX; | dT OT / AX,)° 


Equations (8) and (10) give explicit expressions for p, U; and p in terms of the arbitrary 
vector A; . A special case may be obtained by choosing A; = —0¢/dX; ; this assumes 
[from (8)] that curl pV; = 0, a condition which is always true if the flow depends only 
upon a single space variable and the time, but otherwise is restrictive. Then 


a” : ‘ 
U; = -#¢_/¥ ?, a -V’"d, y= ~<8 + ®, (11) 
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where ® is determined in terms of ¢ by 


o® te) { Oo 0 2 | ' 

ax, ~ ax, | OX aT’ OX; ST /v i (12) 
Expressions (11) and (12) are McVittie’s results but with the dependence of p on $ given 
explicitly. It is necessary to point out, however, that neither Eq. (10) nor Eq. (12) 
has a solution unless the right-hand side is an irrotational vector field. Since (12) is 
equivalent to MeVittie’s equations, this shows that these too have solutions only for a 
restricted class of ¢. 

It is now shown how (3), (4), (5), (6), (7) may be reduced to (11) and (12). To elimi- 
nate ¢, , de ¢; in (5), (6), (7), an expression for the term 4 b & 1 (Vo; — 0°¢;/dX;) in 
(5) is required in terms of ¢. However, if the three quantities white h are equated in (7) 
are denoted by EF, , E£, and £;, respectively, and 4(/, + EF, + E;) is denoted by Y, it 
it seen that (5) isp = —0°¢/dT” + VW; it remains, therefore, to demonstrate that ¥ = ®. 
Differentiating E, with respect to Bs we have 


F an , a2 ‘ ( 42 2 ) 
ov 0k 0 do ; 0 3 g j og ace ‘ 
kee eo eS ae ee oo ae eee ee oe Vor. (13) 
aX OX 0X; OX, OX, OX 0X, \\dX, 07 ) 
Expressions for 0°¢,/0X,0X, and 0°o;/9X.0X, are given in terms of @ by (6); hence 
ov /aX, can be given in terms of ¢ alone. We have 


ny ‘) j Ad a” = ) re) f 7° a | 
— == a. a —_ = | vs 
OX, OX~s (OAs OT OX. OT 7 OX. \d0X_ OT ‘ax, oT 


Fs) S/ Od 
ax, Be oT )/ v6} oF 


a j ra} d Ob 2 \ 
aX, \aX,; dT aX | ¥ ut (14) 


av /aX, and dV/dX; are obtained similarly and we observe that the expressions agree 


with (12); hence, ® and VW are indeed the same. 


BOOK REVIEWS 


An introduction to the theory of diffe rential equations. By Walter Leighton. McGraw-Hill 
Book Company, Inc., New York, Toronto, London, 1952. viii + 174 pp. $3.50 


This book constitutes a carefully written exposition of the elements of the theory of ordinary linear 
differential equations. Accordingly, considerable attention has been given to questions of rigor, and to the 
understanding of basic concepts as opposed to mere formal facility of operation. To this end the author 


discusses the nature of solutions of differential equations and states and discusses existence theorems at 
appropriate points in the text. The proofs of these and other of the more difficult theorems are given, but 
1 appenaices tor the sake of cl arity of exposition. 

The reviewer feels, however, that this book will be too difficult for the student who has just emerged 
from the first course in the calculus as given in most of our colleges. Frequent use is made of results ordi- 
narily not met until a second calculus course is taken, although, where this is done, a statement to this 
effect is made or a reference given to material in the appendices. In several instances results are borrowed 
from more advanced differential equation theory, and in /two cases reference is made to the theory of 
differential equations in the complex domain. This is presumably done for the purpose of presenting a 
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less adumbrated picture of the theory than that usually given, but does not seem likely to contribute 
much to the student’s understanding at this stage. 

With the exceptions noted below, the topics treated are much the same as those in most elementary 
texts on differential equations. The manner in which they are treated, however, is more satisfying to the 
mathematically minded than is usually the case, and the omission of some of the special types of equations 
and the formalism for treating them makes for a more coherent book. The applications are principally 
relegated to a chapter on particle mechanics. This would seem to indicate that this text is aimed at the 
mathematics rather than the engineering student, and although the author does not say so, the general 
tenor of the book bears this out. The better engineering student will nevertheless find much of profit to 
him from the point of view adopted in this book. 

Topics treated which are considered cursorily or not at all in most first texts in differential equations 
are: the second order equation with a regular singular point, the method of successive approximations, 
eigenvalues and eigenfunctions, expansion of functions in series of orthogonal functions, and, in particular, 
oscillation theorems. The treatment of eigenfunctions, eigenvalues, and the related problem of expansion 
of functions in series of orthogonal functions is, to the reviewer, too brief. The discussion on oscillation 
theorems is a distinct novelty in a book on this level. 

The main purpose served by the referencing of advanced material and the appendices would seem 
to be that of stimulating the superior student (for whom the book is most suitable). Would it not then 
be suitable to mention the Lipschitz condition in connection with the existence theorem in Chapter 1? 
Similarly, would it not be stimulating to bring in the notion of vectors and matrices in the Chapter on 
systems of linear differential equations? 

This text is one of the publisher’s International Series in Pure and Applied Mathematics, and is a 
worthy addition thereto. 

WiiuraMm H. Pei 


An introduction to relaxation methods. By F.S. Shaw. Dover Publications, Inc., New York, 
1953. 396 pp. $5.50. 


This book presents an easily read account of how to solve problems by the Relaxation Method. 
Although many papers and several books have been published during the 20 or so years that have elapsed 
since R. V. Southwell first exploited the Relaxation Method of solution none have succeeded as well as 
the present book in presenting to the reader so clear an account of what to do and how to do it. 

A brief historical introduction is followed by chapters dealing with successively more difficult prob- 
lems. A discussion of the solution of Linear Algebraic Equations is followed in succession by a treatment 
of Linear Ordinary Differential Equations, Laplace’s and Poisson’s equations for rectangular boundaries 
and curved boundaries, other more general second order partial differential equations, and fourth order 
partial differential equations. In Chapter VIII the author presents in detail the method of solving Eigen- 
value problems, selecting illustrations from algebraic systems, and ordinary, and partial differential 
equations. In the last chapter (IX) a brief discussion is presented of problems involving an unknown 
boundary (a plastic torsion problem given as illustration), the precision of relaxation solutions and how 
to improve it on a given net, and some dimensional considerations. This last chapter covers too much too 
rapidly but is probably the best compromise between too long a book and no mention of these topics. 
It is perhaps permissible to hope that by the time the reader arrives at the last chapter a briefer treatment 
will serve to extend his knowledge into these additional areas. 

The author disavows any attempt to discuss the physical significance (except in the historical intro- 
duction) of the equation being solved and of the method of solution itself. As a consequence an occasional 
discussion is more general and abstract than necessary. Unfortunately several such items appear in the 
second chapter in the discussion of finite difference approximations. Any reader finding some of the earlier 
sections difficult can safely skip over them to the chapters discussing specific problems since the latter are 
reasonably self contained and the problems solved are treated in detail. It is the reviewer's belief that a 
moderate amount of discussion of physical significance would have increased the readability and under- 
standability of the book. 

Anyone wishing to learn the relaxation method in order to solve problems in this way when appropri- 
ate would do well to study this book and complete the examples used by the author. 

Howarp W. Emmons 
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Foundations of the nonlinear theory of elasticity. By V. V. Novozhilov. Translated from 
the first (1948) Russian edition by F. Bagemihl, H. Komm and W. Seidel. Graylock 
Press, Rochester, N. Y. 1953. vi + 233 pp. $4.00. 


The author’s object is to present an “. . . exposition of the theory of elasticity without any assump- 
tions restricting the magnitude of elongations, displacements or angles of rotation.”’ The first five chapters 
are concerned with the development of the kinematics of finite deformations and the derivation of ‘‘stress- 
strain” relations, equations of motion and boundary conditions. The remaining two chapters discuss the 
problem of elastic stability and the deformation of flexible bodies from the standpoint of the theory 
developed in the earlier chapters. Much of the content of these two chapters is new to the reviewer and 
they form significant contributions to the subjects they discuss. 

Throughout the book, the formulae and equations are expressed in a notation which would have been 
more acceptable to Cauchy than it will be to the modern American reader. 

Taken as a whole, the book cannot be regarded as a comprehensive account of the present state of 
development of finite elasticity theory. This cannot be held entirely to the account of the author, since 
a number of advances have been made after 1948, the year of publication of the original Russian edition. 
These have been critically discussed by C. Truesdell (J. Rat’l Mechanics and Analysis, 1, 125-300 (1952) ; 
2, 593-616 (1953)). 


R. S. Riviirn 


Stochastic processes. By J. L. Doob, John Wiley & Sons, Inc., New York, and Chapman 
& Hall, Ltd., London, 1953. vii + 654 pp. $10.00. 


This book contains the first complete and detailed treatment of stochastic processes and is essential 
for anyone concerned with the mathematical basis of the subject. 

The first two chapters deal with the necessary probability background and definition of the various 
classes of stochastic processes—Gaussian, Markov, stationary, martingales, etc. There follow detailed 
sections on processes with mutually independent, uncorrelated and orthogonal random variables, discrete 
and continuous parameter Markov processes, martingales, processes with independent and orthogonal 
increments and discrete and continuous parameter stationary processes. The final chapter treats the 
theory of linear least squares prediction. 

There is a supplement at the end of the book in which the author outlines those aspects of measure 
theory most needed for the subject matter of the book. However, a fairly thorough knowledge of measure 
theory may nevertheless be considered to be an essential prerequisite for the profitable use of this work. 

LeonaRD C. MaxrImMon 


Theory of matrices. By Sam Perlis. Addison-Wesley Press, Inc., Cambridge 42, Mass., 
1952. xiv + 237 pp. $5.50. 


The book gives a concise and clear treatment of the theory of matrices with emphasis on the basic 
ideas rather than particular applications. It can serve well as a text for a year course, especially if supple- 
mented by applications of the material presented in the book. The approach of the text emphasizes 
strongly the algebra of matrices, with chapters on vector spaces, equivalence, rank and inverses, congru- 
ences and Hermitian congruences. Other topics considered include determinants, polynomials over a 
field and matrices with polynomial elements, similarity, characteristic roots and vectors, and linear 
transformations. Exercises are given at the end of each chapter. 

The book provides a good background for anyone who is to specialize in either pure mathematics or 
physics and applied mathematics. Although some prior acquaintance with the basic notions of modern 
algebra would help the student using this text, the concepts are not presented in a form so abstract that 


the student without such background will feel lost. 
LEONARD C. MAxIMON 




















